truncls3.m

书籍“Regularization tools for training large feed-forward neural networks using Automatic Differentiat

function [ x, istop, itn, nrmgrad] = ...
         truncLS3(VEC, m, n, aprodname, mu, D, b, damp, atol, btol, conlim, ...
	 itnlim,w1,b1,a1,deriv1,df1,w2,b2,a2,deriv2,df2,w3,b3,a3,deriv3,df3,P,T)
%        [ x, istop, itn, nrmgrad ] = ...
%        truncLS3(VEC, m, n, aprodname, mu, D, b, damp, atol, btol, conlim, itnlim,...
%	 w1,b1,a1,deriv1,w2,b2,a2,deriv2,w3,b3,a3,deriv3,P,T)
%
%truncLS3 uses a truncation technique to approximately solve
% (   A)D*D(inv)x = b   in the least squares sence if damp = 0 and m>=n.
% (mu*I)
% where D is diagonal scaling matrix.
% D holds the diagonal elements.

% 		Ax = b  if m = n and 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,mu,D),
% 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 C. C. Paige and M. A. Saunders:
% 1. LSQR: An algorithm for sparse linear equations and sparse least squares,
%    ACM TOMS 8, 1 (Mar 1982) 43-71.
% 2. Algorithm 583.  LSQR: Sparse linear equations and least squares problems,
%    ACM TOMS 8, 2 (Jun 1982) 195-209.
%
% Other potential output parameters:
%    anorm, acond, rnorm, 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.

%     Initialize.

msg=['The exact solution is  x = 0                              '
     'Ax - b is small enough, given atol, btol                  '
     'The least-squares solution is good enough, given atol     '
     'The estimate of cond(Abar) has exceeded conlim            '
     'Ax - b is small enough for this machine                   '
     'The least-squares solution is good enough for this machine'
     'Cond(Abar) seems to be too large for this machine         '
     'The iteration limit has been reached                      '];

itn    = 0;		istop  = 0;		nstop  = 0;
ctol   = 0;		if conlim > 0, ctol = 1/conlim; end;
anorm  = 0;		acond  = 0;
dampsq = damp^2;	ddnorm = 0;		res2   = 0;
xnorm  = 0;		xxnorm = 0;		z      = 0;
cs2    = -1;		sn2    = 0;
major_it=VEC(1);
termind=VEC(2);


% Set up the first vectors u and v for the bidiagonalization.
% These satisfy  beta*u = b,  alfa*v = A'u.

u      = b(1:m);	x    = zeros(n,1);
alfa   = 0;		beta = norm( u );
if beta > 0

   u = (1/beta) * u;	v = D .* feval( aprodname, 2, m, n, u, mu, D, ...
   w1,b1,a1,deriv1,df1,w2,b2,a2,deriv2,df2,w3,b3,a3,deriv3,df3,P,T);
   alfa = norm( v );
%**************
nrmgrad=alfa*beta;
eta=min(1/major_it, nrmgrad);
%nrmrhs=alfa;
end
if alfa > 0
   v = (1/alfa) * v;    w = v;
end

arnorm = alfa * beta;   if arnorm == 0, disp(msg(1,:)); return, end

rhobar = alfa;		phibar = beta;		bnorm  = beta;
rnorm  = beta;
head1  = '   Itn      x(1)      Function';
head2  = ' Compatible   LS      Norm A   Cond A';

%disp([head1 head2])
test1  = 1;		test2  = alfa / beta;
str1   = sprintf( '%6g %12.5e %10.3e',   itn, x(1), rnorm );
str2   = sprintf( '  %8.1e %8.1e',       test1, test2 );
%disp([str1 str2])

%------------------------------------------------------------------
%     Main iteration loop.
%------------------------------------------------------------------
while itn < itnlim
      itn = itn + 1;
%     Perform the next step of the bidiagonalization to obtain the
%     next  beta, u, alfa, v.  These satisfy the relations
%                beta*u  =  a*v   -  alfa*u,
%                alfa*v  =  A'*u  -  beta*v.

%***************
      u    = feval( aprodname, 1, m, n, D .* v, mu, D, ...
      w1,b1,a1,deriv1,df1,w2,b2,a2,deriv2,df2,w3,b3,a3,deriv3,df3,P,T)  -  alfa*u;
      beta = norm( u );
      if beta > 0
         u     = (1/beta) * u;
         anorm = norm([anorm alfa beta damp]);
%****************
         v     = D .* feval( aprodname, 2, m, n, u, mu, D, ...
	 w1,b1,a1,deriv1,df1,w2,b2,a2,deriv2,df2,w3,b3,a3,deriv3,df3,P,T)  -  beta*v;
         alfa  = norm( v );
         if alfa > 0,  v = (1/alfa) * v; end
      end

%     Use a plane rotation to eliminate the damping parameter.
%     This alters the diagonal (rhobar) of the lower-bidiagonal matrix.

      rhobar1 = norm([rhobar damp]);
      cs1     = rhobar / rhobar1;
      sn1     = damp   / rhobar1;
      psi     = sn1 * phibar;
      phibar  = cs1 * phibar;

%     Use a plane rotation to eliminate the subdiagonal element (beta)
%     of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.

      rho     = norm([rhobar1 beta]);
      cs      =   rhobar1/ rho;
      sn      =   beta   / rho;
      theta   =   sn * alfa;
      rhobar  = - cs * alfa;
      phi     =   cs * phibar;
      phibar  =   sn * phibar;
      tau     =   sn * phi;

%     Update x and w.

      ddnorm  = ddnorm + (( 1 / rho)*norm( w ))^2;
      x       = x      + (phi / rho)*w;
      w       = v      - (theta/rho)*w;

%     Use a plane rotation on the right to eliminate the
%     super-diagonal element (theta) of the upper-bidiagonal matrix.
%     Then use the result to estimate  norm(x).

      delta   =   sn2 * rho;
      gambar  = - cs2 * rho;
      rhs     =   phi  -  delta * z;
      zbar    =   rhs / gambar;
      xnorm   =   sqrt(xxnorm + zbar^2);
      gamma   =   norm([gambar theta]);
      cs2     =   gambar / gamma;
      sn2     =   theta  / gamma;
      z       =   rhs    / gamma;
      xxnorm  =   xxnorm  +  z^2;

%     Test for convergence.
%     First, estimate the condition of the matrix  Abar,
%     and the norms of  rbar  and  Abar'rbar.

      acond   =   anorm * sqrt( ddnorm );
      res1    =   phibar^2;
      res2    =   res2  +  psi^2;
      rnorm   =   sqrt( res1 + res2 );
      arnorm  =   alfa * abs( tau );

%     Now use these norms to estimate certain other quantities,
%     some of which will be small near a solution.

      test1   =   rnorm / bnorm;
RNORM(itn,1)=rnorm;
RNORM(itn,2)=test1;
RNORM(itn,3)=xnorm;
RNORM(itn,4)=arnorm;
RNORM(itn,5)=arnorm/bnorm;
RNORM(itn,6)=arnorm/nrmgrad;
RNORM(itn,7)=eta;
      test2   =   arnorm/( anorm * rnorm );
      test3   =       1 / acond;
      t1      =   test1 / (1    +  anorm * xnorm / bnorm);
      rtol    =   btol  +  atol *  anorm * xnorm / bnorm;

%     The following tests guard against extremely small values of
%     atol, btol  or  ctol.  (The user may have set any or all of
%     the parameters  atol, btol, conlim  to 0.)
%     The effect is equivalent to the normal tests using
%     atol = eps,  btol = eps,  conlim = 1/eps.

%if (RNORM(itn,5) < 1e-3) & (RNORM(itn,3) < 1e-2),   istop = 8; end
%if RNORM(itn,5) < 0.5, istop = 8; end

if termind == 1
	if arnorm < nrmgrad*eta, istop = 9; end
	if itn >= itnlim,   istop = 7; end
elseif termind == 3
	if itn >= itnlim,   istop = 7; end
else %termind == 2
	if itn >= itnlim,   istop = 7; end
      if 1 + test3  <= 1, istop = 6; end
      if 1 + test2  <= 1, istop = 5; end
      if 1 + t1     <= 1, istop = 4; end

%     Allow for tolerances set by the user.

      if  test3 <= ctol,  istop = 3; end
      if  test2 <= atol,  istop = 2; end
      if  test1 <= rtol,  istop = 1; end
end
if istop > 0, break, end

%     See if it is time to print something.

      prnt = 0;
      if n     <= 40       , prnt = 1; end
      if itn   <= 10       , prnt = 1; end
      if itn   >= itnlim-10, prnt = 1; end
      if rem(itn,10) == 0  , prnt = 1; end
      if test3 <=  2*ctol  , prnt = 1; end
      if test2 <= 10*atol  , prnt = 1; end
      if test1 <= 10*rtol  , prnt = 1; end
      if istop ~=  0       , prnt = 1; end

      if prnt == 1
         str1 = sprintf( '%6g %12.5e %10.3e',   itn, x(1), rnorm );
         str2 = sprintf( '  %8.1e %8.1e',       test1, test2 );
         str3 = sprintf( ' %8.1e %8.1e',        anorm, acond );
         %disp([str1 str2 str3])
      end
      if istop > 0, break, end
end
x=D .* x;
%[istop itn bnorm nrmgrad]
%RNORM
%pause

%     End of iteration loop.
%     Print the stopping condition.

%disp(' ')
%disp('LSQR finished')
%disp(msg(istop+1,:))
%disp(' ')
str1 = sprintf( 'istop  =%8g   itn    =%8g',      istop, itn    );
str2 = sprintf( 'anorm  =%8.1e   acond  =%8.1e',  anorm, acond  );
str3 = sprintf( 'rnorm  =%8.1e   arnorm =%8.1e',  rnorm, arnorm );
str4 = sprintf( 'bnorm  =%8.1e   xnorm  =%8.1e',  bnorm, xnorm  );
%disp([str1 '   ' str2])
%disp([str3 '   ' str4])

%     End of LSQR