libgenua
Basic Geometry, Numerical Algorithms and Interfaces

LSMR: Iterative solver for leastsquares 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.
X = LSMR(A,B) solves the system of linear equations A*X=B. If the system is inconsistent, it solves the leastsquares problem min b  Ax_2. A is a rectangular matrix of dimension mbyn, 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 leastsquares 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.0e6 (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.0e6 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 = 1e7. 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 leastsquares 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 illconditioned 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 (vvectors generated by the GolubKahan 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 leastsquares 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(BA*X).
[X,ISTOP,ITN,NORMR,NORMAR] = LSMR(A,B,...) gives an estimate of the residual for the normal equation: NORMAR = NORM(A'*(BA*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 matrixvector products A*v and A'*u. For further information, see D. C.L. Fong and M. A. Saunders, LSMR: An iterative algorithm for sparse leastsquares 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 Chinlung Fong clfon Institute for Computational and Mathematical Engineering Stanford University g@st anfor d.ed u
Michael Saunders saund Systems Optimization Laboratory Dept of MS&E, Stanford University. ers@ stanf ord. edu
#include <lsmr.h>
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 = 1e6 
convergence criteria  
double  m_btol = 1e6 
convergence criteria  
double  m_conlim = 1e8 
condition number limit (actually, limit for the estimated condition)  
size_t  m_maxiter = 128 
permitted number of iteration  