Revision 8cd68f3f650d00762f17885aeb455c7b4f8c0ac2 authored by Michael Friedlander on 19 September 2014, 17:06:45 UTC, committed by Michael Friedlander on 19 September 2014, 17:06:45 UTC
lsqr.m
function [ x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var ]...
= lsqr( m, n, A, b, damp, atol, btol, conlim, itnlim, show )
%
% [ x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var ]...
% = lsqr( m, n, A, b, damp, atol, btol, conlim, itnlim, show );
%
% 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 or a function handle of aprod( mode,x ),
% 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.
%-----------------------------------------------------------------------
% 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:
% 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.
%
% 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.
% r1norm = norm(r), where r = b - Ax.
% r2norm = sqrt( norm(r)^2 + damp^2 * norm(x)^2 )
% = r1norm if damp = 0.
% anorm = estimate of Frobenius norm of Abar = [ A ].
% [damp*I]
% acond = estimate of cond(Abar).
% arnorm = estimate of norm(A'*r - damp^2*x).
% xnorm = norm(x).
% var (if present) estimates all diagonals of (A'A)^{-1} (if damp=0)
% or more generally (A'A + damp^2*I)^{-1}.
% This is well defined if A has full column rank or damp > 0.
% (Not sure what var means if rank(A) < n and damp = 0.)
%
%
% 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.
% 07 Aug 2002: Output parameter rnorm replaced by r1norm and r2norm.
% Michael Saunders, Systems Optimization Laboratory,
% Dept of MS&E, Stanford University.
% 03 Jul 2007: Modified 'aprodname' to A, which can either be an m by n
% matrix, or a function handle.
% Ewout van den Berg, University of British Columbia
% 03 Jul 2007: Modified 'test2' condition, omitted 'test1'.
% Ewout van den Berg, University of British Columbia
%-----------------------------------------------------------------------
% 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 '];
wantvar= nargout >= 6;
if wantvar, var = zeros(n,1); end
if show
disp(' ')
disp('LSQR Least-squares solution of Ax = b')
str1 = sprintf('The matrix A has %8g rows and %8g cols', m, n);
str2 = sprintf('damp = %20.14e wantvar = %8g', damp,wantvar);
str3 = sprintf('atol = %8.2e conlim = %8.2e', atol, conlim);
str4 = sprintf('btol = %8.2e itnlim = %8g' , btol, itnlim);
disp(str1); disp(str2); disp(str3); disp(str4);
end
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;
% 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 = Aprod(u,2);
alfa = norm( v );
end
if alfa > 0
v = (1/alfa) * v; w = v;
end
arnorm = alfa * beta;
if arnorm == 0
if show, disp(msg(1,:)); end
return
end
arnorm0= arnorm;
rhobar = alfa; phibar = beta; bnorm = beta;
rnorm = beta;
r1norm = rnorm;
r2norm = rnorm;
head1 = ' Itn x(1) r1norm r2norm ';
head2 = ' Compatible LS Norm A Cond A';
if show
disp(' ')
disp([head1 head2])
test1 = 1; test2 = alfa / beta;
str1 = sprintf( '%6g %12.5e', itn, x(1) );
str2 = sprintf( ' %10.3e %10.3e', r1norm, r2norm );
str3 = sprintf( ' %8.1e %8.1e', test1, test2 );
disp([str1 str2 str3])
end
%------------------------------------------------------------------
% 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 = Aprod(v,1) - alfa*u;
beta = norm( u );
if beta > 0
u = (1/beta) * u;
anorm = norm([anorm alfa beta damp]);
v = Aprod(u, 2) - 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.
t1 = phi /rho;
t2 = - theta/rho;
dk = (1/rho)*w;
x = x + t1*w;
w = v + t2*w;
ddnorm = ddnorm + norm(dk)^2;
if wantvar, var = var + dk.*dk; end
% 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 );
% 07 Aug 2002:
% Distinguish between
% r1norm = ||b - Ax|| and
% r2norm = rnorm in current code
% = sqrt(r1norm^2 + damp^2*||x||^2).
% Estimate r1norm from
% r1norm = sqrt(r2norm^2 - damp^2*||x||^2).
% Although there is cancellation, it might be accurate enough.
r1sq = rnorm^2 - dampsq * xxnorm;
r1norm = sqrt( abs(r1sq) ); if r1sq < 0, r1norm = - r1norm; end
r2norm = rnorm;
% Now use these norms to estimate certain other quantities,
% some of which will be small near a solution.
test1 = rnorm / bnorm;
test2 = arnorm / arnorm0;
% 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 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
% 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
if show
str1 = sprintf( '%6g %12.5e', itn, x(1) );
str2 = sprintf( ' %10.3e %10.3e', r1norm, r2norm );
str3 = sprintf( ' %8.1e %8.1e', test1, test2 );
str4 = sprintf( ' %8.1e %8.1e', anorm, acond );
disp([str1 str2 str3 str4])
end
end
if istop > 0, break, end
end
% End of iteration loop.
% Print the stopping condition.
if show
disp(' ')
disp('LSQR finished')
disp(msg(istop+1,:))
disp(' ')
str1 = sprintf( 'istop =%8g r1norm =%8.1e', istop, r1norm );
str2 = sprintf( 'anorm =%8.1e arnorm =%8.1e', anorm, arnorm );
str3 = sprintf( 'itn =%8g r2norm =%8.1e', itn, r2norm );
str4 = sprintf( 'acond =%8.1e xnorm =%8.1e', acond, xnorm );
disp([str1 ' ' str2])
disp([str3 ' ' str4])
disp(' ')
end
%-----------------------------------------------------------------------
% End of lsqr.m
%-----------------------------------------------------------------------
function z = Aprod(x,mode)
if mode == 1
if isnumeric(A), z = A*x;
else z = A(x,1);
end
else
if isnumeric(A), z = (x'*A)';
else z = A(x,2);
end
end
end % function Aprod
end
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...