libgenua Basic Geometry, Numerical Algorithms and Interfaces
LsmrSolver< Scalar > Class Template Reference

## Detailed Description

### template<typename Scalar> class LsmrSolver< Scalar >

LSMR: Iterative solver for least-squares problems.

The operator 'A' passed to solve() needs to provide the interface

```   void muladd(const DVector<Scalar> &x, DVector<Scalar> &y) const;
void muladdTransposed(const DVector<Scalar> &x, DVector<Scalar> &y) const;```

which compute y = A*x and y = transpose(A)*x.

CraigSolver

X = LSMR(A,B) solves the system of linear equations A*X=B. If the system is inconsistent, it solves the least-squares problem min ||b - Ax||_2. A is a rectangular matrix of dimension m-by-n, where all cases are allowed: m=n, m>n, or m<n. B is a vector of length m. The matrix A may be dense or sparse (usually sparse).

X = LSMR(AFUN,B) takes a function handle AFUN instead of the matrix A. AFUN(X,1) takes a vector X and returns A*X. AFUN(X,2) returns A'*X. AFUN can be used in all the following syntaxes.

X = LSMR(A,B,LAMBDA) solves the regularized least-squares problem min ||(B) - ( A )X|| ||(0) (LAMBDA*I) ||_2 where LAMBDA is a scalar. If LAMBDA is [] or 0, the system is solved without regularization.

X = LSMR(A,B,LAMBDA,ATOL,BTOL) continues iterations until a certain backward error estimate is smaller than some quantity depending on ATOL and BTOL. Let RES = B - A*X be the residual vector for the current approximate solution X. If A*X = B seems to be consistent, LSMR terminates when NORM(RES) <= ATOL*NORM(A)*NORM(X) + BTOL*NORM(B). Otherwise, LSMR terminates when NORM(A'*RES) <= ATOL*NORM(A)*NORM(RES). If both tolerances are 1.0e-6 (say), the final NORM(RES) should be accurate to about 6 digits. (The final X will usually have fewer correct digits, depending on cond(A) and the size of LAMBDA.) If ATOL or BTOL is [], a default value of 1.0e-6 will be used. Ideally, they should be estimates of the relative error in the entries of A and B respectively. For example, if the entries of A have 7 correct digits, set ATOL = 1e-7. This prevents the algorithm from doing unnecessary work beyond the uncertainty of the input data.

X = LSMR(A,B,LAMBDA,ATOL,BTOL,CONLIM) 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. If CONLIM is [], the default value is CONLIM = 1e+8. Maximum precision can be obtained by setting ATOL = BTOL = CONLIM = 0, but the number of iterations may then be excessive.

X = LSMR(A,B,LAMBDA,ATOL,BTOL,CONLIM,ITNLIM) terminates if the number of iterations reaches ITNLIM. The default is ITNLIM = min(m,n). For ill-conditioned systems, a larger value of ITNLIM may be needed.

X = LSMR(A,B,LAMBDA,ATOL,BTOL,CONLIM,ITNLIM,LOCALSIZE) runs LSMR with reorthogonalization on the last LOCALSIZE v_k's (v-vectors generated by the Golub-Kahan bidiagonalization). LOCALSIZE = 0 or [] runs LSMR without reorthogonalization. LOCALSIZE = Inf specifies full reorthogonalization of the v_k's. Reorthogonalizing only u_k or both u_k and v_k are not an option here, because reorthogonalizing all v_k's makes the u_k's close to orthogonal. Details are given in the submitted SIAM paper.

X = LSMR(A,B,LAMBDA,ATOL,BTOL,CONLIM,ITNLIM,LOCALSIZE,SHOW) prints an iteration log if SHOW=true. The default value is SHOW=false.

[X,ISTOP] = LSMR(A,B,...) gives the reason for termination. ISTOP = 0 means X=0 is a solution. = 1 means X is an approximate solution to A*X = B, according to ATOL and BTOL. = 2 means X approximately solves the least-squares problem according to ATOL. = 3 means COND(A) seems to be greater than CONLIM. = 4 is the same as 1 with ATOL = BTOL = EPS. = 5 is the same as 2 with ATOL = EPS. = 6 is the same as 3 with CONLIM = 1/EPS. = 7 means ITN reached ITNLIM before the other stopping conditions were satisfied.

[X,ISTOP,ITN] = LSMR(A,B,...) gives ITN = the number of LSMR iterations.

[X,ISTOP,ITN,NORMR] = LSMR(A,B,...) gives an estimate of the residual norm: NORMR = norm(B-A*X).

[X,ISTOP,ITN,NORMR,NORMAR] = LSMR(A,B,...) gives an estimate of the residual for the normal equation: NORMAR = NORM(A'*(B-A*X)).

[X,ISTOP,ITN,NORMR,NORMAR,NORMA] = LSMR(A,B,...) gives an estimate of the Frobenius norm of A.

[X,ISTOP,ITN,NORMR,NORMAR,NORMA,CONDA] = LSMR(A,B,...) gives an estimate of the condition number of A.

[X,ISTOP,ITN,NORMR,NORMAR,NORMA,CONDA,NORMX] = LSMR(A,B,...) gives an estimate of NORM(X).

LSMR uses an iterative method requiring matrix-vector products A*v and A'*u. For further information, see D. C.-L. Fong and M. A. Saunders, LSMR: An iterative algorithm for sparse least-squares problems, SIAM J. Sci. Comput., submitted 1 June 2010. See http://www.stanford.edu/~clfong/lsmr.html.

08 Dec 2009: First release version of LSMR. 09 Apr 2010: Updated documentation and default parameters. 14 Apr 2010: Updated documentation. 03 Jun 2010: LSMR with local and/or full reorthogonalization of v_k. 10 Mar 2011: Bug fix in reorthgonalization. (suggested by David Gleich)

David Chin-lung Fong clfon.nosp@m.g@st.nosp@m.anfor.nosp@m.d.ed.nosp@m.u Institute for Computational and Mathematical Engineering Stanford University

Michael Saunders saund.nosp@m.ers@.nosp@m.stanf.nosp@m.ord..nosp@m.edu Systems Optimization Laboratory Dept of MS&E, Stanford University.

`#include <lsmr.h>`

Inheritance diagram for LsmrSolver< Scalar >:
[legend]
Collaboration diagram for LsmrSolver< Scalar >:
[legend]

## Public Member Functions

void reortho (size_t n)
set the number of reorthogonalization vectors

template<class Operator >
ExitCode solve (const Operator &A, const DVector< Scalar > &b, DVector< Scalar > &x, Scalar lambda=0)
minimize |Ax - b|^2 + lambda |x|^2

Public Member Functions inherited from SolIterativeSolver
void tolerance (int maxiter, double atoler, double btoler, double conlim=1e8)
set convergence tolerances and maximum number of iterations

## Private Attributes

size_t m_vcols = 16
maximum dimension of subspace V

## Additional Inherited Members

Static Public Member Functions inherited from SolIterativeSolver
static const char * statusMessage (int code)
access text corresponding to error code

static bool success (int code)
does exit code indicate succes (in some sense, at least)?

Protected Attributes inherited from SolIterativeSolver
double m_atol = 1e-6
convergence criteria

double m_btol = 1e-6
convergence criteria

double m_conlim = 1e8
condition number limit (actually, limit for the estimated condition)

size_t m_maxiter = 128
permitted number of iteration

The documentation for this class was generated from the following file: