LSQR method (for linear systems and linear least-squares problems). More...

#include <BelosLSQRSolMgr.hpp>

Inheritance diagram for Belos::LSQRSolMgr< ScalarType, MV, OP, scalarTypeIsComplex >:

Public Member Functions
	LSQRSolMgr ()

	LSQRSolMgr (const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)

virtual	~LSQRSolMgr ()

Teuchos::RCP< SolverManager < ScalarType, MV, OP > >	clone () const override
	clone for Inverted Injection (DII) More...

Private Types
typedef Details::RealSolverManager < ScalarType, MV, OP, isComplex >	base_type

Static Private Attributes
static const bool	isComplex = Teuchos::ScalarTraits<ScalarType>::isComplex

Detailed Description

template<class ScalarType, class MV, class OP, const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex>
class Belos::LSQRSolMgr< ScalarType, MV, OP, scalarTypeIsComplex >

LSQR method (for linear systems and linear least-squares problems).

Author: Sarah Knepper and David Day

Template Parameters

ScalarType	The type of entries in the right-hand side vector(s) and solution vector(s) .
MV	The multivector type; the type of the solution vector(s) and right-hand side vector(s).
OP	The type of the matrix (and any preconditioner, if one is provided).

Algorithm

LSQR (Paige and Saunders; see References) is an iterative method for solving linear least-squares problems and linear systems. It can solve any of the following problems:

Solve for
Find that minimizes $\|Ax - b\|_2^2$
Find that minimizes $\|Ax - b\|_2^2 + \lambda^2 \|x\|_2^2$

The third problem above is the most general and includes the previous two. This is the problem LSQR actually solves. Here, $\lambda$ is a user-provided positive real constant (the "damping parameter") which regularizes the problem so that it always has a bounded solution, even if $A$ does not have full rank.

In the words of Paige and Saunders: "The method is based on the Golub-Kahan bidiagonalization process. It is algebraically equivalent to applying MINRES to the normal equation[s] $(A^T A + \lambda 2I)x = A^T b$ , but has better numerical properties, especially if $A$ is ill-conditioned."

LSQR has some special algorithmic properties:

It reduces $\|b - A x\|_2$ (the two-norm of the residual) monotonically.
LSQR also computes a monotonically increasing estimate of the two-norm condition number of the matrix .

Property #2 makes LSQR useful for mixed-precision algorithms. If the matrix $A$ has condition number greater than the inverse of machine precision in the current working precision, one can reconstruct the problem to solve in the next higher precision and restart, possibly using the previous solution as an initial guess.

ScalarType must be real

This LSQR implementation currently only supports real-valued (not complex-valued) ScalarType types. You may check whether ScalarType is complex using the following code:

if (Teuchos::ScalarTraits<ScalarType>::isComplex) {
 // ScalarType is complex valued.
} else {
 // ScalarType is real valued.
}

This is not a limitation of the LSQR method itself, just of the current implementation. If there is sufficient interest, we can remedy this deficiency. For now, if you attempt to invoke the constructor when ScalarType is complex, the constructor will throw an exception. This is why this class inherits from Details::RealSolverManager. LSQRSolMgr can still compile if ScalarType is complex, but you will not be able to construct a LSQRSolMgr instance in that case, due to the aforementioned run-time error that the constructor raises. We do this so that the class will still compile, whether ScalarType is real or complex. This helps make SolverFactory valid to compile, whether ScalarType is real or complex.

Preconditioning

If the linear problem to solve includes a preconditioner (in the LinearProblem object), then the least-squares problem is solved for the preconditioned linear system. Preconditioning changes the least-squares problem (in the sense of changing the norms), and the solution depends on the preconditioner in this sense. In the context of linear least-squares problems, "preconditioning" refers to the regularization matrix. In this solver, the regularization matrix is always a scalar multiple of the identity (standard form least squares).

A converged preconditioned residual norm suffices for convergence, but is not necessary. LSQR sometimes returns a larger relative residual norm than what would have been returned by a linear solver. For details on the stopping criteria, see the documentation of LSQRStatusTest, which implements the three-part stopping criterion recommended by Paige and Saunders.

Some Belos solvers implement detection of "loss of accuracy." That refers to the difference between convergence of the original linear system and convergence of the (left-)preconditioned linear system. LSQR does not implement detection of "loss of accuracy," because it is unclear what this means for linear least squares in general. This LSQR solves a possibly inconsistent system in a least-squares sense.

References

C. C. Paige and M. A. Saunders, LSQR: An algorithm for sparse linear equations and sparse least squares, TOMS 8(1), 43-71 (1982).

C. C. Paige and M. A. Saunders, Algorithm 583; LSQR: Sparse linear equations and least-squares problems, TOMS 8(2), 195-209 (1982).