Lsqr

aprod(mode, m, n, x, iw, rw)

This is the simplest example for testing LSQR. \(A = rw\). If mode = 1, aprod computes \(y = A x\). Ff mode = 2, aprod computes \(y = A^T x\). for some matrix A.

cgls(Aname, shift, b, m, n, kmax, tol, prnt)

Solves a symmetric system

\((A^T A + shift I)x = A^T b\) or \(N x = A^T b\),

where A is a general m by n matrix and shift may be positive or negative. The method should be stable if \(N = (A^T A + shift I)\) is positive definite. It MAY be unstable otherwise.

‘gemat’ specifies an M-file gemat.m (say), such that:

  • y = gemat(0, x, m, n) identifies gemat (if you wish),

  • y = gemat(1, x, m, n) gives \(y = A x\),

  • y = gemat(2, x, m, n) gives \(y = A^T x\).

The M-file must not be called Aname.m!

USAGE:

[x, resNE, k, info] = cgls( Aname ,shift, b, m, n, kmax, tol, prnt)

INPUTS:

Aname: name of file shift: if 0 then cgls is Hestenes and Stiefel’s specialized form b: used in the formulas in description m, n: dimensions of the matrix kmax: maximum number of iterations. tol: desired relative residual size, \(norm(rNE)/norm(A^T b)\),

where \(rNE = A^T b - N x\).

prnt: 1 gives an iteration log, 0 suppresses it.

OUTPUTS:

x: from the formulas in description resNE: relative residual for normal equations: \(norm(A^T b - Nx)/norm(A^T b)\). k: final number of iterations performed. info: gives information on the result:

  • 1 if required accuracy was achieved,

  • 2 if the limit kmax was reached,

  • 3 if N seems to be singular or indefinite.

  • 4 if instability seems likely. (N indefinite and normx decreased.)

Note

When \(shift = 0\), cgls is Hestenes and Stiefel’s specialized form of the conjugate-gradient method for least-squares problems. The general shift is a simple modification.

lsqr(m, n, aprodname, b, damp, atol, btol, conlim, itnlim, show, covindex, userstop)

LSQR solves \(A x = b\) or \(min ||b - A x||_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), where aprod is a function handle that performs the matrix-vector operations. If mode = 1, aprod must return \(y = A x\) without altering x. If mode = 2, aprod must return \(y = A^T x\) without altering x.

USAGE:

[x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var, cov, userout] = lsqr( m, n, aprodname, b, damp, atol, btol, conlim, itnlim, show, covindex, userstop)

INPUTS:

m,n: dimensions of the matrix aprodname: function handle b: from \(A x = b\) 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, 0 - suppresses output, covindex: specified index set of the columns and rows of the

covariance matrix

userstop: a function handle that checks a user-defined stopping

criteria. The function handle should be in the form [istop, atol, userout] = userstop(istop, atol, arnorm, itn);

if userstop returns istop > 0 lsqr will stop if userstop returns istop = 0 lsqr will continue the userstop function can alter the tolerance atol.

OUTPUTS:

x: is the final solution. istop: gives the reason for termination.

  • 1 means x is an approximate solution to Ax = b.

  • 2 means x approximately solves the least-squares problem.

itn: changed to 0 in the code r1norm: = norm(r), where \(r = b - A x\). 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^T r - damp^2 x)\). xnorm: = norm(x). var: (if present) estimates all diagonals of \((A^T A)^{-1}\) (if \(damp = 0\))

or more generally \((A^T 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\).)

cov: (if present) esitmates the rows and columns of the

covariance matrix specified by the covindex set.

userout: output data from the userstop function

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.

lsqraprod(mode, m, n, x, iw, rw)

if mode = 1, lsqraprod computes \(y = A x\), if mode = 2, lsqraprod computes \(y = A^T x\) for some matrix A.

This is a simple example for testing LSQR. It uses the leading m n submatrix from:

Example

A = [ 1
1 2
2 3
3 4
…

n ]

suitably padded by zeros.

lsqrtest(m, n, damp)

If \(m = n\) and \(damp = 0\), this sets up a system \(A x = b\) and calls lsqr.m to solve it. Otherwise, the usual least-squares or damped least-squares problem is solved. lsqraprod.m defines the m x n matrix A.