pdsco.m

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

%-----------------------------------------------------------------------
%  Backtracking linesearch.
%-----------------------------------------------------------------------
   fail     =  true;
   nf       =  0;

   while nf < maxf
      nf       =   nf + 1;
      xnew     =   x  +  stepx * dx;
      ynew     =   y  +  stepz * dy;
      znew     =   z  +  stepz * dz;
      [obj,grad,hess] = feval( Fname, (xnew*beta) );
      obj      =   obj /theta;
      grad     =   grad*(beta /theta)  +  gamma2*xnew;
      Q        =   hess*(beta2/theta)  +  gamma2;
      rlinnew  =   rlin     -  stepx * Adx;
      tlinnew  =   tlin     -  stepz *(Ady + dz);
      rnew     =   rlinnew  -  LSdamp2 * ynew;
      tnew     =   tlinnew  +  grad;
      Xznew    =   xnew .* znew;
      vnew     =   mu  -  Xznew;
      fnew     =   norm([rnew;  tnew;  vnew]);     % Must be 2-norm here.
      step     =   min( stepx, stepz );
      % if fnew >= f,
      %    fratio = fnew / f;
      %    fprintf( 'fnew / f%8.1e', fratio)
      % end
      if fnew <= (1 - eta*step)*f
         fail = false;
         break;
      end

%     If the first attempt doesn't work,
%     make stepx and stepz the same.

      if nf == 1 & stepx ~= stepz
         stepx   = min( stepx, stepz );
         stepz   = stepx;
      elseif nf < maxf
         stepx   = stepx/2;
         stepz   = stepx;
      end;
   end

   if fail
      fprintf('\n     Linesearch failed (nf too big)');
      nfail = nfail + 1;
   else
      nfail = 0;
   end

   x = xnew;   r = rnew;   rlin = rlinnew;   f  = fnew;
   y = ynew;   t = tnew;   tlin = tlinnew;   Xz = Xznew;
   z = znew;   v = vnew;

%-----------------------------------------------------------------------
%  Set convergence measures.
%-----------------------------------------------------------------------

   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 );
   primalfeas    = Pinf  <=  featol;
   dualfeas      = Dinf  <=  featol;
   complementary = Cinf  <=  opttol;
   enough        = PDitns>=       4;  % Prevent premature termination.
   converged     = primalfeas  &  dualfeas  &  complementary  &  enough;

%-----------------------------------------------------------------------
%  Iteration log.
%-----------------------------------------------------------------------
   str1    = sprintf('\n%3g%5.1f' , PDitns     , log10(mu)   );
   str2    = sprintf('%6.3f%6.3f' , stepx      , stepz       );
   if stepx < 0.001 | stepz < 0.001
      str2 = sprintf('%6.1e%6.1e' , stepx      , stepz       );
   end

   str3    = sprintf('%6.1f%6.1f' , log10(Pinf), log10(Dinf));
   str4    = sprintf('%6.1f%15.7e', log10(Cinf), objtrue     );
   str5    = sprintf('%3g%8.1f'   , nf         , center      );
   if center > 99999
      str5 = sprintf('%3g%8.1e'   , nf         , center      );
   end
   fprintf([str1 str2 str3 str4 str5])
   if direct
      % relax
   else
      fprintf(' %5.1f%7g%7.3f', log10(atolold), itncg, r3ratio)
   end


%-----------------------------------------------------------------------
%  Test for termination.
%-----------------------------------------------------------------------
   if kminor
      fprintf( 'Start of next minor itn...\n')
      keyboard
   end

   if converged
      fprintf('\n   Converged')
   elseif PDitns >= maxitn
      fprintf('\n   Too many iterations')
      inform = 1;
      break
   elseif nfail  >= maxfail
      fprintf('\n   Too many linesearch failures')
      inform = 2;
      break
   else

      % Reduce mu and LSdamp, and reset certain residuals.

      stepmu  = min( stepx , stepz   );
      stepmu  = min( stepmu, steptol );
      muold   = mu;
      mu      = mu   -  stepmu * mu;
      if center >= 100,  mu = muold;  end

%     mutrad  = mu0 * (sum(Xz) / n); % 24 May 1998:  Traditional value. But--
%     mu      = min( mu, mutrad );   % it seemed to decrease mu too much.

      mu      = max( mu, mulast );   % 13 Jun 1998:  No need for smaller mu.
      LSdamp  = max( sqrt(mu), delta );
      LSdamp  = min( LSdamp  , delta0);
      LSdamp2 = LSdamp^2;
      r       = rlin  -  LSdamp2 * y;
      v       = mu  -  Xz;
      f       = norm([r; t; v]);

%     Reduce atol for LSQR (and SYMMLQ).
%     NOW DONE AT TOP OF LOOP.

      atolold = atol;
      % if atol > atol2
      %   atolfac = (mu/mufirst)^0.25;
      %   atol    = max( atol*atolfac, atol2 );
      % end

%     atol    = min( atol, mu );     % 22 Jan 2001:  a la Inexact Newton.
%     atol    = min( atol, 0.5*mu ); % 30 Jan 2001:  A bit tighter

      % If the linesearch took more than one function (nf > 1),
      % we assume the search direction needed more accuracy
      % (though this may be true only for LPs).
      % 12 Jun 1998: Ask for more accuracy if nf > 2.
      % 24 Nov 2000: Also if the steps are small.
      % 30 Jan 2001: Small steps might be ok with warm start.
      % 06 Feb 2001: Not necessarily.  Reinstated tests in next line.

      if nf > 2  |  min( stepx, stepz ) <= 0.01
         atol = atolold*0.1;
      end
   end
end
%-----------------------------------------------------------------------
% End of main loop.
%-----------------------------------------------------------------------

fprintf('\n\nmax( x ) =%10.3f', max(x)     )
fprintf('    max(|y|) =%10.3f', max(abs(y)))
fprintf('    max( z ) =%10.3f', max(z)     )
fprintf('   scaled')

x = x*beta;   y = y*zeta;   z = z*zeta;   % Unscale x, y, z.

fprintf(  '\nmax( x ) =%10.3f', max(x)     )
fprintf('    max(|y|) =%10.3f', max(abs(y)))
fprintf('    max( z ) =%10.3f', max(z)     )
fprintf(' unscaled')

time   = cputime - time;
str1   = sprintf('\nPDitns   =%10g', PDitns );
str2   = sprintf(     ' itns =%10g', CGitns );
fprintf( [str1 ' ' solver str2] )
fprintf('    time     =%10.1f', time);

pdsxxxdistrib(x,z);   % Private function

if wait
   keyboard
end
%-----------------------------------------------------------------------
% End function pdsco.m
%-----------------------------------------------------------------------


function pdsxxxdistrib(x,z)

% pdsxxxdistrib(x) or pdsxxxdistrib(x,z) prints the
% distribution of 1 or 2 vectors.
%
% 18 Dec 2000.  First version with 2 vectors.

two  = nargin > 1;
fprintf('\n\nDistribution of vector     x')
if two, fprintf('         z'); end

x1   = 10^(floor(log10(max(x))) + 1);
z1   = 10^(floor(log10(max(z))) + 1);
x1   = max(x1,z1);
kmax = 10;

for k = 1:kmax
    x2 = x1;    x1 = x1/10;
    if k==kmax, x1 = 0; end
    nx = length(find(x>=x1 & x<x2));
    fprintf('\n[%7.3g,%7.3g )%10g', x1, x2, nx);
    if two
       nz = length(find(z>=x1 & z<x2));
       fprintf('%10g', nz);
    end
end

disp(' ')

%-----------------------------------------------------------------------
% End private function pdsxxxdistrib
%-----------------------------------------------------------------------


function [ x, istop, itn, outfo ] = ...
   pdsxxxlsqr( m, n, aprodname, iw, rw, b, damp, ...
	       atol, btol, conlim, itnlim, show, info )

% Special version of LSQR for use with pdsco.m.
% It continues with a reduced atol if a pdsco-specific test isn't
% satisfied with the input atol.
%
% LSQR solves  Ax = b  or  min ||b - Ax||_2  if damp = 0,
% or   min || (b)  -  (  A   )x ||   otherwise.
%          || (0)     (damp I)  ||2
% A  is an m by n matrix defined by  y = aprod( mode,m,n,x,iw,rw ),
% where the parameter 'aprodname' refers to a function 'aprod' that
% performs the matrix-vector operations.
% If mode = 1,   aprod  must return  y = Ax   without altering x.
% If mode = 2,   aprod  must return  y = A'x  without altering x.
% WARNING:   The file containing the function 'aprod'
%            must not be called aprodname.m !!!!

%-----------------------------------------------------------------------
% LSQR uses an iterative (conjugate-gradient-like) method.
% For further information, see 
% 1. C. C. Paige and M. A. Saunders (1982a).
%    LSQR: An algorithm for sparse linear equations and sparse least squares,
%    ACM TOMS 8(1), 43-71.
% 2. C. C. Paige and M. A. Saunders (1982b).
%    Algorithm 583.  LSQR: Sparse linear equations and least squares problems,
%    ACM TOMS 8(2), 195-209.
% 3. M. A. Saunders (1995).  Solution of sparse rectangular systems using
%    LSQR and CRAIG, BIT 35, 588-604.
%
% Input parameters:
% iw, rw      are not used by lsqr, but are passed to aprod.
% atol, btol  are stopping tolerances.  If both are 1.0e-9 (say),
%             the final residual norm should be accurate to about 9 digits.
%             (The final x will usually have fewer correct digits,
%             depending on cond(A) and the size of damp.)
% conlim      is also a stopping tolerance.  lsqr terminates if an estimate
%             of cond(A) exceeds conlim.  For compatible systems Ax = b,
%             conlim could be as large as 1.0e+12 (say).  For least-squares
%             problems, conlim should be less than 1.0e+8.
%             Maximum precision can be obtained by setting
%             atol = btol = conlim = zero, but the number of iterations
%             may then be excessive.
% itnlim      is an explicit limit on iterations (for safety).
% show = 1    gives an iteration log,
% show = 0    suppresses output.
% info        is a structure special to pdsco.m, used to test if
%             was small enough, and continuing if necessary with smaller atol.
%
%
% Output parameters:
% x           is the final solution.
% istop       gives the reason for termination.
% istop       = 1 means x is an approximate solution to Ax = b.
%             = 2 means x approximately solves the least-squares problem.
% rnorm       = norm(r) if damp = 0, where r = b - Ax,
%             = sqrt( norm(r)**2  +  damp**2 * norm(x)**2 ) otherwise.
% xnorm       = norm(x).
% var         estimates diag( inv(A'A) ).  Omitted in this special version.
% outfo       is a structure special to pdsco.m, returning information
%             about whether atol had to be reduced.
%             
% Other potential output parameters:
% anorm, acond, arnorm, xnorm
%
%        1990: Derived from Fortran 77 version of LSQR.
% 22 May 1992: bbnorm was used incorrectly.  Replaced by anorm.
% 26 Oct 1992: More input and output parameters added.
% 01 Sep 1994: Matrix-vector routine is now a parameter 'aprodname'.
%              Print log reformatted.
% 14 Jun 1997: show  added to allow printing or not.
% 30 Jun 1997: var   added as an optional output parameter.
%              It returns an estimate of diag( inv(A'A) ).
% 12 Feb 2001: atol  can now be reduced and iterations continued if necessary.
%              info, outfo are new problem-dependent parameters for such purposes.
%              In this version they are specialized for pdsco.m.
%-----------------------------------------------------------------------

%     Initialize.

msg=['The exact solution is  x = 0                              '
     'Ax - b is small enough, given atol, btol                  '