Eigen/src/IterativeLinearSolvers/LSMR.h - mirror - Git at Google


            
            // This file is part of Eigen, a lightweight C++ template library

            
            // for linear algebra.

            
            //

            
            // This Source Code Form is subject to the terms of the Mozilla

            
            // Public License v. 2.0. If a copy of the MPL was not distributed

            
            // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

            
            // SPDX-FileCopyrightText: The Eigen Authors

            
            // SPDX-License-Identifier: MPL-2.0


            #ifndef EIGEN_LSMR_H

            
            #define EIGEN_LSMR_H


            // LSMR is an iterative algorithm for least-squares problems  min ||A x - b||,

            
            // based on the Golub-Kahan bidiagonalization process.  It is analytically

            
            // equivalent to applying MINRES to the normal equation  A^T A x = A^T b, so the

            
            // quantity ||A^T r_k|| decreases monotonically (where r_k = b - A x_k).  In

            
            // practice ||r_k|| also decreases monotonically, which makes LSMR safer than

            
            // LSQR to stop early.  With a damping parameter lambda > 0 it instead minimizes

            
            // || (A; lambda I) x - (b; 0) ||, i.e. it solves the regularized (Tikhonov)

            
            // least-squares problem.

            
            //

            
            // This implementation follows the published algorithm:

            
            //

            
            //   D. C.-L. Fong and M. A. Saunders, "LSMR: An Iterative Algorithm for Sparse

            
            //   Least-Squares Problems", SIAM J. Sci. Comput. 33(5):2950-2971, 2011.

            
            //   https://web.stanford.edu/group/SOL/software/lsmr/

            
            //

            
            // The scalar recurrences, the cheap estimates of ||r||, ||A^T r||, ||A|| and

            
            // cond(A), and the stopping rules reproduce the reference Fortran 90

            
            // implementation by the same authors (distributed under the BSD / Common Public

            
            // License at the URL above).  No source code from that implementation is

            
            // reproduced here; only the mathematical algorithm is implemented, in Eigen's

            
            // own idiom.


            // IWYU pragma: private

            
            #include "./InternalHeaderCheck.h"


            namespace Eigen {


            namespace internal {


            /** \internal Low-level LSMR algorithm

            
             * \param mat       The matrix A.

            
             * \param rhs       The right hand side vector b.

            
             * \param x         On input the initial guess x0 (usually zero), on output the

            
             *                  computed solution.

            
             * \param precond   A preconditioner.  It is applied as a self-adjoint right

            
             *                  preconditioner: LSMR is run on the better-conditioned system

            
             *                  \f$ A M^{-1} z = b \f$ (where \c precond.solve() applies

            
             *                  \f$ M^{-1} \f$) and the solution is recovered as

            
             *                  \f$ x = M^{-1} z \f$.  With IdentityPreconditioner this

            
             *                  reduces exactly to the reference algorithm.

            
             * \param iters     On input the maximum number of iterations, on output the

            
             *                  number of iterations performed.

            
             * \param tol_error On output, an estimate of the relative residual of the

            
             *                  normal equations \f$ ||A^T r|| / (||A||\,||r||) \f$.

            
             * \param atol      Stopping tolerance bounding the assumed relative error in

            
             *                  the entries of \a A.

            
             * \param btol      Stopping tolerance bounding the assumed relative error in

            
             *                  the entries of \a b.

            
             * \param lambda    The damping parameter \f$ \lambda \ge 0 \f$ (0 for an

            
             *                  unregularized problem).

            
             * \param conlim    An upper limit on cond(A); iterations stop if the estimate

            
             *                  exceeds it.  Pass 0 to disable (equivalent to 1/eps).

            
             *

            
             * \returns the LSMR stopping code \c istop:

            
             *   - 0: x = 0 (i.e. the supplied guess) is the exact solution; no iterations.

            
             *   - 1: A x = b is compatible and ||r|| is small enough, given atol and btol.

            
             *   - 2: a least-squares solution is good enough, given atol.

            
             *   - 3: the estimate of cond(A) exceeded conlim.

            
             *   - 4: A x = b is compatible and ||r|| is as small as machine precision allows.

            
             *   - 5: a least-squares solution is as good as machine precision allows.

            
             *   - 6: cond(A) seems too large for this machine.

            
             *   - 7: the iteration limit was reached.

            
             */

            
            template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>

            
            EIGEN_DONT_INLINE Index lsmr(const MatrixType& mat, const Rhs& rhs, Dest& x, const Preconditioner& precond,

            
                                         Index& iters, typename Dest::RealScalar& tol_error, const typename Dest::RealScalar& atol,

            
                                         const typename Dest::RealScalar& btol, const typename Dest::RealScalar& lambda,

            
                                         const typename Dest::RealScalar& conlim) {

            
              using numext::abs;

            
              using numext::sqrt;

            
              using RealScalar = typename Dest::RealScalar;

            
              using Scalar = typename Dest::Scalar;

            
              using VectorType = Matrix<Scalar, Dynamic, 1>;


              const RealScalar zero(0);

            
              const RealScalar one(1);


              const Index n = mat.cols();

            
              const Index maxIters = iters;


              // n-vectors needed before the early-return below. u is the only m-vector.

            
              VectorType v(n), Atu(n);


              // Set up the first vectors u and v for the bidiagonalization. These satisfy

            
              // beta*u = b - A*x0  and  alpha*v = M^{-1} A^T u.

            
              VectorType u = rhs - mat * x;  // working residual r0 of the initial guess

            
              RealScalar alpha = zero;

            
              RealScalar beta = u.stableNorm();

            
              if (beta > zero) {

            
                u /= beta;

            
                Atu.noalias() = mat.adjoint() * u;

            
                v = precond.solve(Atu);

            
                alpha = v.stableNorm();

            
                if (alpha > zero) v /= alpha;

            
              }


              iters = 0;

            
              // If b - A*x0 = 0 or A^T(b - A*x0) = 0 then the current x already solves the

            
              // (least-squares) problem: no correction is needed.

            
              if (alpha * beta == zero) {

            
                tol_error = zero;

            
                return 0;  // istop = 0

            
              }


              // dx accumulates the correction (in z-space when a preconditioner is used).

            
              // The remaining n-vectors are only needed once we start iterating.

            
              VectorType dx = VectorType::Zero(n);

            
              VectorType h = v;

            
              VectorType hbar = VectorType::Zero(n);

            
              VectorType t(n);


              // Quantities driving the two plane rotations.

            
              RealScalar alphabar = alpha;

            
              RealScalar zetabar = alpha * beta;

            
              RealScalar rho = one;

            
              RealScalar rhobar = one;

            
              RealScalar cbar = one;

            
              RealScalar sbar = zero;


              // Quantities for the running estimate of ||r||.

            
              RealScalar betadd = beta;

            
              RealScalar betad = zero;

            
              RealScalar rhodold = one;

            
              RealScalar tautildeold = zero;

            
              RealScalar thetatilde = zero;

            
              RealScalar zeta = zero;

            
              RealScalar d = zero;


              // Quantities for the running estimates of ||A|| and cond(A).

            
              RealScalar normA2 = alpha * alpha;

            
              RealScalar maxrbar = zero;

            
              RealScalar minrbar = NumTraits<RealScalar>::highest();


              const RealScalar normb = beta;

            
              const RealScalar ctol = conlim > zero ? one / conlim : zero;

            
              RealScalar test2 = zero;  // recomputed each iteration; also read after the loop for tol_error


              Index istop = 0;

            
              Index itn = 0;

            
              while (istop == 0) {

            
                ++itn;


                // Perform the next step of the bidiagonalization to obtain the next

            
                // beta, u, alpha, v.  These satisfy

            
                //     beta*u = A*(M^{-1} v) - alpha*u,

            
                //     alpha*v = M^{-1} A^T u - beta*v.

            
                t = precond.solve(v);

            
                u *= -alpha;

            
                u.noalias() += mat * t;

            
                beta = u.stableNorm();

            
                if (beta > zero) {

            
                  u /= beta;

            
                  Atu.noalias() = mat.adjoint() * u;

            
                  t = precond.solve(Atu);

            
                  v = t - beta * v;

            
                  alpha = v.stableNorm();

            
                  if (alpha > zero) v /= alpha;

            
                }


                // Construct rotation Qhat_{k} that folds in the damping.

            
                const RealScalar alphahat = numext::hypot(alphabar, lambda);

            
                const RealScalar chat = alphabar / alphahat;

            
                const RealScalar shat = lambda / alphahat;


                // Use a plane rotation Q_{k} to turn B_{k} into R_{k}.

            
                const RealScalar rhoold = rho;

            
                rho = numext::hypot(alphahat, beta);

            
                const RealScalar c = alphahat / rho;

            
                const RealScalar s = beta / rho;

            
                const RealScalar thetanew = s * alpha;

            
                alphabar = c * alpha;


                // Use a plane rotation Qbar_{k} to turn R_{k}^T into Rbar_{k}.

            
                const RealScalar rhobarold = rhobar;

            
                const RealScalar zetaold = zeta;

            
                const RealScalar thetabar = sbar * rho;

            
                const RealScalar rhotemp = cbar * rho;

            
                rhobar = numext::hypot(cbar * rho, thetanew);

            
                cbar = cbar * rho / rhobar;

            
                sbar = thetanew / rhobar;

            
                zeta = cbar * zetabar;

            
                zetabar = -sbar * zetabar;


                // Update h, hbar and the (correction) solution dx.

            
                hbar = h - (thetabar * rho / (rhoold * rhobarold)) * hbar;

            
                dx += (zeta / (rho * rhobar)) * hbar;

            
                h = v - (thetanew / rho) * h;


                // Estimate ||r||.  Apply rotation Qhat_{k}, then Q_{k}, then Qtilde_{k-1}.

            
                const RealScalar betaacute = chat * betadd;

            
                const RealScalar betacheck = -shat * betadd;

            
                const RealScalar betahat = c * betaacute;

            
                betadd = -s * betaacute;


                const RealScalar thetatildeold = thetatilde;

            
                const RealScalar rhotildeold = numext::hypot(rhodold, thetabar);

            
                const RealScalar ctildeold = rhodold / rhotildeold;

            
                const RealScalar stildeold = thetabar / rhotildeold;

            
                thetatilde = stildeold * rhobar;

            
                rhodold = ctildeold * rhobar;

            
                betad = -stildeold * betad + ctildeold * betahat;


                tautildeold = (zetaold - thetatildeold * tautildeold) / rhotildeold;

            
                const RealScalar taud = (zeta - thetatilde * tautildeold) / rhodold;

            
                d += betacheck * betacheck;

            
                const RealScalar normr = sqrt(d + numext::abs2(betad - taud) + numext::abs2(betadd));


                // Estimate ||A||.

            
                normA2 += beta * beta;

            
                const RealScalar normA = sqrt(normA2);

            
                normA2 += alpha * alpha;


                // Estimate cond(A).

            
                maxrbar = numext::maxi(maxrbar, rhobarold);

            
                if (itn > 1) minrbar = numext::mini(minrbar, rhobarold);

            
                const RealScalar condA = numext::maxi(maxrbar, rhotemp) / numext::mini(minrbar, rhotemp);


                // Compute the norms needed for the stopping rules.

            
                const RealScalar normAr = abs(zetabar);

            
                const RealScalar normx = dx.stableNorm();


                const RealScalar test1 = normr / normb;

            
                test2 = (normA * normr > zero) ? normAr / (normA * normr) : zero;

            
                const RealScalar test3 = one / condA;

            
                const RealScalar t1 = test1 / (one + normA * normx / normb);

            
                const RealScalar rtol = btol + atol * normA * normx / normb;


                // The "1 + test <= 1" guards trigger near machine precision and make the

            
                // method behave as if atol = btol = eps and conlim = 1/eps even when the

            
                // user passed 0 for any of them.  The user tolerances are tested first so a

            
                // genuine convergence (istop 1/2/3) takes priority over the machine limits.

            
                // (istop 4 uses t1 rather than test1, matching the reference algorithm.)

            
                if (test1 <= rtol)

            
                  istop = 1;

            
                else if (test2 <= atol)

            
                  istop = 2;

            
                else if (test3 <= ctol)

            
                  istop = 3;

            
                else if (one + t1 <= one)

            
                  istop = 4;

            
                else if (one + test2 <= one)

            
                  istop = 5;

            
                else if (one + test3 <= one)

            
                  istop = 6;

            
                else if (itn >= maxIters)

            
                  istop = 7;

            
              }


              // Recover the solution: x <- x0 + M^{-1} dx.

            
              t = precond.solve(dx);

            
              x += t;


              iters = itn;

            
              tol_error = test2;

            
              return istop;

            
            }


            }  // namespace internal


            template <typename MatrixType_, typename Preconditioner_ = IdentityPreconditioner>

            
            class LSMR;


            namespace internal {


            template <typename MatrixType_, typename Preconditioner_>

            
            struct traits<LSMR<MatrixType_, Preconditioner_> > {

            
              using MatrixType = MatrixType_;

            
              using Preconditioner = Preconditioner_;

            
            };


            }  // namespace internal


            /** \ingroup IterativeLinearSolvers_Module

            
             * \brief An LSMR solver for sparse (or dense) least-squares problems

            
             *

            
             * This class solves for the least-squares solution of \c A \c x = \c b using the

            
             * LSMR algorithm of Fong and Saunders.  LSMR is based on the Golub-Kahan

            
             * bidiagonalization and is analytically equivalent to MINRES applied to the

            
             * normal equation \f$ A^T A x = A^T b \f$, so the residual of the normal

            
             * equation \f$ ||A^T r|| \f$ decreases monotonically.  In exact arithmetic LSMR

            
             * returns the minimum-norm least-squares solution.  The matrix \c A can be

            
             * non-symmetric and rectangular; \c A and the vectors \c x and \c b can be

            
             * either dense or sparse.  LSMR only needs \c A through the products \f$ Av \f$

            
             * and \f$ A^T u \f$, so a matrix-free operator may be used as well.

            
             *

            
             * Unlike LeastSquaresConjugateGradient, which forms \f$ A^T A \f$ implicitly,

            
             * LSMR works on \c A directly through the bidiagonalization and is therefore

            
             * more robust on ill-conditioned problems.

            
             *

            
             * \tparam MatrixType_ the type of the matrix A, can be a dense or a sparse matrix.

            
             * \tparam Preconditioner_ the type of the preconditioner. Default is IdentityPreconditioner.

            
             *

            
             * \implsparsesolverconcept

            
             *

            
             * The maximum number of iterations and the tolerance can be controlled via the

            
             * setMaxIterations() and setTolerance() methods. The defaults are twice the

            
             * number of columns of the matrix for the maximum number of iterations and

            
             * NumTraits<Scalar>::epsilon() for the tolerance. setTolerance() sets both of

            
             * the algorithm's stopping tolerances \c atol (relative error assumed in \c A)

            
             * and \c btol (relative error assumed in \c b); they can also be set

            
             * independently via setToleranceA() and setToleranceB().

            
             *

            
             * Unlike most other iterative solvers, error() does not report the relative

            
             * residual \f$ ||Ax-b||/||b|| \f$: it reports the estimate

            
             * \f$ ||A^T r|| / (||A||\,||r||) \f$, with \f$ r = b - Ax \f$, of the relative

            
             * residual of the normal equations, the quantity that the least-squares

            
             * stopping rule bounds by \c atol.

            
             *

            
             * The setDamping() method enables Tikhonov regularization: with a damping

            
             * \f$ \lambda > 0 \f$ the solver minimizes

            
             * \f$ ||Ax-b||^2 + \lambda^2 ||x||^2 \f$, for which a unique solution always

            
             * exists. The setConditionLimit() method can be used to stop the iterations as

            
             * soon as the estimated condition number of \c A exceeds a given bound.

            
             *

            
             * This class can be used like the other iterative solvers. Here is a typical

            
             * usage example:

            
               \code

            
               int m = 1000000, n = 10000;

            
               VectorXd x(n), b(m);

            
               SparseMatrix<double> A(m, n);

            
               // fill A and b

            
               LSMR<SparseMatrix<double> > lsmr;

            
               lsmr.compute(A);

            
               x = lsmr.solve(b);

            
               std::cout << "#iterations:     " << lsmr.iterations() << std::endl;

            
               std::cout << "estimated error: " << lsmr.error()      << std::endl;

            
               // update b, and solve again

            
               x = lsmr.solve(b);

            
               \endcode

            
             *

            
             * By default the iterations start with x=0 as an initial guess of the solution.

            
             * One can control the start using the solveWithGuess() method.

            
             *

            
             * If a non-default preconditioner is supplied it is applied as a self-adjoint

            
             * right preconditioner (LSMR is run on \f$ A M^{-1} z = b \f$ and the solution

            
             * recovered from \f$ M x = z \f$); it should therefore be self-adjoint, as

            
             * DiagonalPreconditioner and IdentityPreconditioner are.

            
             *

            
             * \sa class LeastSquaresConjugateGradient, class ConjugateGradient, SparseLU, SparseQR

            
             */

            
            template <typename MatrixType_, typename Preconditioner_>

            
            class LSMR : public IterativeSolverBase<LSMR<MatrixType_, Preconditioner_> > {

            
             protected:

            
              using Base = IterativeSolverBase<LSMR>;

            
              using Base::m_error;

            
              using Base::m_info;

            
              using Base::m_isInitialized;

            
              using Base::m_iterations;

            
              using Base::matrix;


             public:

            
              using MatrixType = MatrixType_;

            
              using Scalar = typename MatrixType::Scalar;

            
              using RealScalar = typename MatrixType::RealScalar;

            
              using Preconditioner = Preconditioner_;


              /** Default constructor. */

            
              LSMR() : Base() {}


              /** Initialize the solver with matrix \a A for further \c Ax=b solving.

            
               *

            
               * This constructor is a shortcut for the default constructor followed

            
               * by a call to compute().

            
               *

            
               * \warning this class stores a reference to the matrix A as well as some

            
               * precomputed values that depend on it. Therefore, if \a A is changed

            
               * this class becomes invalid. Call compute() to update it with the new

            
               * matrix A, or modify a copy of A.

            
               */

            
              template <typename MatrixDerived>

            
              explicit LSMR(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}


              /** Sets the damping parameter \f$ \lambda \ge 0 \f$ for Tikhonov

            
               * regularization. With \a lambda > 0 the solver minimizes

            
               * \f$ ||Ax-b||^2 + \lambda^2 ||x||^2 \f$. The default is 0 (no damping).

            
               *

            
               * \note When a non-identity preconditioner \f$ M \f$ is used the damping

            
               * applies in the preconditioned variable, i.e. the solver minimizes

            
               * \f$ ||Ax-b||^2 + \lambda^2 ||Mx||^2 \f$. With the default

            
               * IdentityPreconditioner this is exactly \f$ ||Ax-b||^2 + \lambda^2 ||x||^2 \f$.

            
               */

            
              LSMR& setDamping(const RealScalar& lambda) {

            
                m_lambda = lambda;

            
                return *this;

            
              }


              /** \returns the damping parameter. \sa setDamping() */

            
              RealScalar damping() const { return m_lambda; }


              /** Sets an upper limit on the estimated condition number of \a A. The

            
               * iterations stop as soon as the estimate exceeds \a conlim. The default is

            
               * 0, which disables the limit (equivalent to 1/epsilon).

            
               */

            
              LSMR& setConditionLimit(const RealScalar& conlim) {

            
                m_conditionLimit = conlim;

            
                return *this;

            
              }


              /** \returns the condition-number limit. \sa setConditionLimit() */

            
              RealScalar conditionLimit() const { return m_conditionLimit; }


              /** Sets the stopping tolerance \c atol, which bounds the relative error

            
               * assumed in the entries of \a A. It drives the least-squares stopping rule

            
               * \f$ ||A^T r|| \le atol\,||A||\,||r|| \f$. If left unset (the default) it

            
               * falls back to tolerance(). \sa setToleranceB(), setTolerance() */

            
              LSMR& setToleranceA(const RealScalar& atol) {

            
                m_atol = atol;

            
                return *this;

            
              }


              /** \returns \c atol, or tolerance() if setToleranceA() has not been called.

            
               * \sa setToleranceA() */

            
              RealScalar toleranceA() const { return m_atol >= RealScalar(0) ? m_atol : Base::m_tolerance; }


              /** Sets the stopping tolerance \c btol, which bounds the relative error

            
               * assumed in the entries of \a b. It enters the compatible-system stopping

            
               * rule \f$ ||r|| \le btol\,||b|| + atol\,||A||\,||x|| \f$. If left unset (the

            
               * default) it falls back to tolerance(). \sa setToleranceA(), setTolerance() */

            
              LSMR& setToleranceB(const RealScalar& btol) {

            
                m_btol = btol;

            
                return *this;

            
              }


              /** \returns \c btol, or tolerance() if setToleranceB() has not been called.

            
               * \sa setToleranceB() */

            
              RealScalar toleranceB() const { return m_btol >= RealScalar(0) ? m_btol : Base::m_tolerance; }


              /** \internal */

            
              template <typename Rhs, typename Dest>

            
              void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const {

            
                m_iterations = Base::maxIterations();


                Index istop = internal::lsmr(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error, toleranceA(),

            
                                             toleranceB(), m_lambda, m_conditionLimit);

            
                // istop in {0,1,2,4,5}: the (least-squares) solution was found, possibly

            
                // only to within machine precision (4,5). istop in {3,6,7}: stopped on the

            
                // condition-number limit or the iteration limit without meeting the

            
                // requested tolerance.

            
                m_info = (istop == 3 || istop == 6 || istop == 7) ? NoConvergence : Success;

            
              }


             protected:

            
              RealScalar m_lambda = RealScalar(0);

            
              RealScalar m_conditionLimit = RealScalar(0);

            
              // Negative means "unset": toleranceA()/toleranceB() then fall back to tolerance().

            
              RealScalar m_atol = RealScalar(-1);

            
              RealScalar m_btol = RealScalar(-1);

            
            };


            }  // end namespace Eigen


            #endif  // EIGEN_LSMR_H
	// This file is part of Eigen, a lightweight C++ template library
	// for linear algebra.
	//
	// This Source Code Form is subject to the terms of the Mozilla
	// Public License v. 2.0. If a copy of the MPL was not distributed
	// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
	// SPDX-FileCopyrightText: The Eigen Authors
	// SPDX-License-Identifier: MPL-2.0

	#ifndef EIGEN_LSMR_H
	#define EIGEN_LSMR_H

	// LSMR is an iterative algorithm for least-squares problems min \|\|A x - b\|\|,
	// based on the Golub-Kahan bidiagonalization process. It is analytically
	// equivalent to applying MINRES to the normal equation A^T A x = A^T b, so the
	// quantity \|\|A^T r_k\|\| decreases monotonically (where r_k = b - A x_k). In
	// practice \|\|r_k\|\| also decreases monotonically, which makes LSMR safer than
	// LSQR to stop early. With a damping parameter lambda > 0 it instead minimizes
	// \|\| (A; lambda I) x - (b; 0) \|\|, i.e. it solves the regularized (Tikhonov)
	// least-squares problem.
	//
	// This implementation follows the published algorithm:
	//
	// D. C.-L. Fong and M. A. Saunders, "LSMR: An Iterative Algorithm for Sparse
	// Least-Squares Problems", SIAM J. Sci. Comput. 33(5):2950-2971, 2011.
	// https://web.stanford.edu/group/SOL/software/lsmr/
	//
	// The scalar recurrences, the cheap estimates of \|\|r\|\|, \|\|A^T r\|\|, \|\|A\|\| and
	// cond(A), and the stopping rules reproduce the reference Fortran 90
	// implementation by the same authors (distributed under the BSD / Common Public
	// License at the URL above). No source code from that implementation is
	// reproduced here; only the mathematical algorithm is implemented, in Eigen's
	// own idiom.

	// IWYU pragma: private
	#include "./InternalHeaderCheck.h"

	namespace Eigen {

	namespace internal {

	/** \internal Low-level LSMR algorithm
	* \param mat The matrix A.
	* \param rhs The right hand side vector b.
	* \param x On input the initial guess x0 (usually zero), on output the
	* computed solution.
	* \param precond A preconditioner. It is applied as a self-adjoint right
	* preconditioner: LSMR is run on the better-conditioned system
	* \f$ A M^{-1} z = b \f$ (where \c precond.solve() applies
	* \f$ M^{-1} \f$) and the solution is recovered as
	* \f$ x = M^{-1} z \f$. With IdentityPreconditioner this
	* reduces exactly to the reference algorithm.
	* \param iters On input the maximum number of iterations, on output the
	* number of iterations performed.
	* \param tol_error On output, an estimate of the relative residual of the
	* normal equations \f$ \|\|A^T r\|\| / (\|\|A\|\|\,\|\|r\|\|) \f$.
	* \param atol Stopping tolerance bounding the assumed relative error in
	* the entries of \a A.
	* \param btol Stopping tolerance bounding the assumed relative error in
	* the entries of \a b.
	* \param lambda The damping parameter \f$ \lambda \ge 0 \f$ (0 for an
	* unregularized problem).
	* \param conlim An upper limit on cond(A); iterations stop if the estimate
	* exceeds it. Pass 0 to disable (equivalent to 1/eps).
	*
	* \returns the LSMR stopping code \c istop:
	* - 0: x = 0 (i.e. the supplied guess) is the exact solution; no iterations.
	* - 1: A x = b is compatible and \|\|r\|\| is small enough, given atol and btol.
	* - 2: a least-squares solution is good enough, given atol.
	* - 3: the estimate of cond(A) exceeded conlim.
	* - 4: A x = b is compatible and \|\|r\|\| is as small as machine precision allows.
	* - 5: a least-squares solution is as good as machine precision allows.
	* - 6: cond(A) seems too large for this machine.
	* - 7: the iteration limit was reached.
	*/
	template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
	EIGEN_DONT_INLINE Index lsmr(const MatrixType& mat, const Rhs& rhs, Dest& x, const Preconditioner& precond,
	Index& iters, typename Dest::RealScalar& tol_error, const typename Dest::RealScalar& atol,
	const typename Dest::RealScalar& btol, const typename Dest::RealScalar& lambda,
	const typename Dest::RealScalar& conlim) {
	using numext::abs;
	using numext::sqrt;
	using RealScalar = typename Dest::RealScalar;
	using Scalar = typename Dest::Scalar;
	using VectorType = Matrix<Scalar, Dynamic, 1>;

	const RealScalar zero(0);
	const RealScalar one(1);

	const Index n = mat.cols();
	const Index maxIters = iters;

	// n-vectors needed before the early-return below. u is the only m-vector.
	VectorType v(n), Atu(n);

	// Set up the first vectors u and v for the bidiagonalization. These satisfy
	// betau = b - Ax0 and alpha*v = M^{-1} A^T u.
	VectorType u = rhs - mat * x; // working residual r0 of the initial guess
	RealScalar alpha = zero;
	RealScalar beta = u.stableNorm();
	if (beta > zero) {
	u /= beta;
	Atu.noalias() = mat.adjoint() * u;
	v = precond.solve(Atu);
	alpha = v.stableNorm();
	if (alpha > zero) v /= alpha;
	}

	iters = 0;
	// If b - Ax0 = 0 or A^T(b - Ax0) = 0 then the current x already solves the
	// (least-squares) problem: no correction is needed.
	if (alpha * beta == zero) {
	tol_error = zero;
	return 0; // istop = 0
	}

	// dx accumulates the correction (in z-space when a preconditioner is used).
	// The remaining n-vectors are only needed once we start iterating.
	VectorType dx = VectorType::Zero(n);
	VectorType h = v;
	VectorType hbar = VectorType::Zero(n);
	VectorType t(n);

	// Quantities driving the two plane rotations.
	RealScalar alphabar = alpha;
	RealScalar zetabar = alpha * beta;
	RealScalar rho = one;
	RealScalar rhobar = one;
	RealScalar cbar = one;
	RealScalar sbar = zero;

	// Quantities for the running estimate of \|\|r\|\|.
	RealScalar betadd = beta;
	RealScalar betad = zero;
	RealScalar rhodold = one;
	RealScalar tautildeold = zero;
	RealScalar thetatilde = zero;
	RealScalar zeta = zero;
	RealScalar d = zero;

	// Quantities for the running estimates of \|\|A\|\| and cond(A).
	RealScalar normA2 = alpha * alpha;
	RealScalar maxrbar = zero;
	RealScalar minrbar = NumTraits<RealScalar>::highest();

	const RealScalar normb = beta;
	const RealScalar ctol = conlim > zero ? one / conlim : zero;
	RealScalar test2 = zero; // recomputed each iteration; also read after the loop for tol_error

	Index istop = 0;
	Index itn = 0;
	while (istop == 0) {
	++itn;

	// Perform the next step of the bidiagonalization to obtain the next
	// beta, u, alpha, v. These satisfy
	// betau = A(M^{-1} v) - alpha*u,
	// alphav = M^{-1} A^T u - betav.
	t = precond.solve(v);
	u *= -alpha;
	u.noalias() += mat * t;
	beta = u.stableNorm();
	if (beta > zero) {
	u /= beta;
	Atu.noalias() = mat.adjoint() * u;
	t = precond.solve(Atu);
	v = t - beta * v;
	alpha = v.stableNorm();
	if (alpha > zero) v /= alpha;
	}

	// Construct rotation Qhat_{k} that folds in the damping.
	const RealScalar alphahat = numext::hypot(alphabar, lambda);
	const RealScalar chat = alphabar / alphahat;
	const RealScalar shat = lambda / alphahat;

	// Use a plane rotation Q_{k} to turn B_{k} into R_{k}.
	const RealScalar rhoold = rho;
	rho = numext::hypot(alphahat, beta);
	const RealScalar c = alphahat / rho;
	const RealScalar s = beta / rho;
	const RealScalar thetanew = s * alpha;
	alphabar = c * alpha;

	// Use a plane rotation Qbar_{k} to turn R_{k}^T into Rbar_{k}.
	const RealScalar rhobarold = rhobar;
	const RealScalar zetaold = zeta;
	const RealScalar thetabar = sbar * rho;
	const RealScalar rhotemp = cbar * rho;
	rhobar = numext::hypot(cbar * rho, thetanew);
	cbar = cbar * rho / rhobar;
	sbar = thetanew / rhobar;
	zeta = cbar * zetabar;
	zetabar = -sbar * zetabar;

	// Update h, hbar and the (correction) solution dx.
	hbar = h - (thetabar * rho / (rhoold * rhobarold)) * hbar;
	dx += (zeta / (rho * rhobar)) * hbar;
	h = v - (thetanew / rho) * h;

	// Estimate \|\|r\|\|. Apply rotation Qhat_{k}, then Q_{k}, then Qtilde_{k-1}.
	const RealScalar betaacute = chat * betadd;
	const RealScalar betacheck = -shat * betadd;
	const RealScalar betahat = c * betaacute;
	betadd = -s * betaacute;

	const RealScalar thetatildeold = thetatilde;
	const RealScalar rhotildeold = numext::hypot(rhodold, thetabar);
	const RealScalar ctildeold = rhodold / rhotildeold;
	const RealScalar stildeold = thetabar / rhotildeold;
	thetatilde = stildeold * rhobar;
	rhodold = ctildeold * rhobar;
	betad = -stildeold * betad + ctildeold * betahat;

	tautildeold = (zetaold - thetatildeold * tautildeold) / rhotildeold;
	const RealScalar taud = (zeta - thetatilde * tautildeold) / rhodold;
	d += betacheck * betacheck;
	const RealScalar normr = sqrt(d + numext::abs2(betad - taud) + numext::abs2(betadd));

	// Estimate \|\|A\|\|.
	normA2 += beta * beta;
	const RealScalar normA = sqrt(normA2);
	normA2 += alpha * alpha;

	// Estimate cond(A).
	maxrbar = numext::maxi(maxrbar, rhobarold);
	if (itn > 1) minrbar = numext::mini(minrbar, rhobarold);
	const RealScalar condA = numext::maxi(maxrbar, rhotemp) / numext::mini(minrbar, rhotemp);

	// Compute the norms needed for the stopping rules.
	const RealScalar normAr = abs(zetabar);
	const RealScalar normx = dx.stableNorm();

	const RealScalar test1 = normr / normb;
	test2 = (normA * normr > zero) ? normAr / (normA * normr) : zero;
	const RealScalar test3 = one / condA;
	const RealScalar t1 = test1 / (one + normA * normx / normb);
	const RealScalar rtol = btol + atol * normA * normx / normb;

	// The "1 + test <= 1" guards trigger near machine precision and make the
	// method behave as if atol = btol = eps and conlim = 1/eps even when the
	// user passed 0 for any of them. The user tolerances are tested first so a
	// genuine convergence (istop 1/2/3) takes priority over the machine limits.
	// (istop 4 uses t1 rather than test1, matching the reference algorithm.)
	if (test1 <= rtol)
	istop = 1;
	else if (test2 <= atol)
	istop = 2;
	else if (test3 <= ctol)
	istop = 3;
	else if (one + t1 <= one)
	istop = 4;
	else if (one + test2 <= one)
	istop = 5;
	else if (one + test3 <= one)
	istop = 6;
	else if (itn >= maxIters)
	istop = 7;
	}

	// Recover the solution: x <- x0 + M^{-1} dx.
	t = precond.solve(dx);
	x += t;

	iters = itn;
	tol_error = test2;
	return istop;
	}

	} // namespace internal

	template <typename MatrixType_, typename Preconditioner_ = IdentityPreconditioner>
	class LSMR;

	namespace internal {

	template <typename MatrixType_, typename Preconditioner_>
	struct traits<LSMR<MatrixType_, Preconditioner_> > {
	using MatrixType = MatrixType_;
	using Preconditioner = Preconditioner_;
	};

	} // namespace internal

	/** \ingroup IterativeLinearSolvers_Module
	* \brief An LSMR solver for sparse (or dense) least-squares problems
	*
	* This class solves for the least-squares solution of \c A \c x = \c b using the
	* LSMR algorithm of Fong and Saunders. LSMR is based on the Golub-Kahan
	* bidiagonalization and is analytically equivalent to MINRES applied to the
	* normal equation \f$ A^T A x = A^T b \f$, so the residual of the normal
	* equation \f$ \|\|A^T r\|\| \f$ decreases monotonically. In exact arithmetic LSMR
	* returns the minimum-norm least-squares solution. The matrix \c A can be
	* non-symmetric and rectangular; \c A and the vectors \c x and \c b can be
	* either dense or sparse. LSMR only needs \c A through the products \f$ Av \f$
	* and \f$ A^T u \f$, so a matrix-free operator may be used as well.
	*
	* Unlike LeastSquaresConjugateGradient, which forms \f$ A^T A \f$ implicitly,
	* LSMR works on \c A directly through the bidiagonalization and is therefore
	* more robust on ill-conditioned problems.
	*
	* \tparam MatrixType_ the type of the matrix A, can be a dense or a sparse matrix.
	* \tparam Preconditioner_ the type of the preconditioner. Default is IdentityPreconditioner.
	*
	* \implsparsesolverconcept
	*
	* The maximum number of iterations and the tolerance can be controlled via the
	* setMaxIterations() and setTolerance() methods. The defaults are twice the
	* number of columns of the matrix for the maximum number of iterations and
	* NumTraits<Scalar>::epsilon() for the tolerance. setTolerance() sets both of
	* the algorithm's stopping tolerances \c atol (relative error assumed in \c A)
	* and \c btol (relative error assumed in \c b); they can also be set
	* independently via setToleranceA() and setToleranceB().
	*
	* Unlike most other iterative solvers, error() does not report the relative
	* residual \f$ \|\|Ax-b\|\|/\|\|b\|\| \f$: it reports the estimate
	* \f$ \|\|A^T r\|\| / (\|\|A\|\|\,\|\|r\|\|) \f$, with \f$ r = b - Ax \f$, of the relative
	* residual of the normal equations, the quantity that the least-squares
	* stopping rule bounds by \c atol.
	*
	* The setDamping() method enables Tikhonov regularization: with a damping
	* \f$ \lambda > 0 \f$ the solver minimizes
	* \f$ \|\|Ax-b\|\|^2 + \lambda^2 \|\|x\|\|^2 \f$, for which a unique solution always
	* exists. The setConditionLimit() method can be used to stop the iterations as
	* soon as the estimated condition number of \c A exceeds a given bound.
	*
	* This class can be used like the other iterative solvers. Here is a typical
	* usage example:
	\code
	int m = 1000000, n = 10000;
	VectorXd x(n), b(m);
	SparseMatrix<double> A(m, n);
	// fill A and b
	LSMR<SparseMatrix<double> > lsmr;
	lsmr.compute(A);
	x = lsmr.solve(b);
	std::cout << "#iterations: " << lsmr.iterations() << std::endl;
	std::cout << "estimated error: " << lsmr.error() << std::endl;
	// update b, and solve again
	x = lsmr.solve(b);
	\endcode
	*
	* By default the iterations start with x=0 as an initial guess of the solution.
	* One can control the start using the solveWithGuess() method.
	*
	* If a non-default preconditioner is supplied it is applied as a self-adjoint
	* right preconditioner (LSMR is run on \f$ A M^{-1} z = b \f$ and the solution
	* recovered from \f$ M x = z \f$); it should therefore be self-adjoint, as
	* DiagonalPreconditioner and IdentityPreconditioner are.
	*
	* \sa class LeastSquaresConjugateGradient, class ConjugateGradient, SparseLU, SparseQR
	*/
	template <typename MatrixType_, typename Preconditioner_>
	class LSMR : public IterativeSolverBase<LSMR<MatrixType_, Preconditioner_> > {
	protected:
	using Base = IterativeSolverBase<LSMR>;
	using Base::m_error;
	using Base::m_info;
	using Base::m_isInitialized;
	using Base::m_iterations;
	using Base::matrix;

	public:
	using MatrixType = MatrixType_;
	using Scalar = typename MatrixType::Scalar;
	using RealScalar = typename MatrixType::RealScalar;
	using Preconditioner = Preconditioner_;

	/** Default constructor. */
	LSMR() : Base() {}

	/** Initialize the solver with matrix \a A for further \c Ax=b solving.
	*
	* This constructor is a shortcut for the default constructor followed
	* by a call to compute().
	*
	* \warning this class stores a reference to the matrix A as well as some
	* precomputed values that depend on it. Therefore, if \a A is changed
	* this class becomes invalid. Call compute() to update it with the new
	* matrix A, or modify a copy of A.
	*/
	template <typename MatrixDerived>
	explicit LSMR(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}

	/** Sets the damping parameter \f$ \lambda \ge 0 \f$ for Tikhonov
	* regularization. With \a lambda > 0 the solver minimizes
	* \f$ \|\|Ax-b\|\|^2 + \lambda^2 \|\|x\|\|^2 \f$. The default is 0 (no damping).
	*
	* \note When a non-identity preconditioner \f$ M \f$ is used the damping
	* applies in the preconditioned variable, i.e. the solver minimizes
	* \f$ \|\|Ax-b\|\|^2 + \lambda^2 \|\|Mx\|\|^2 \f$. With the default
	* IdentityPreconditioner this is exactly \f$ \|\|Ax-b\|\|^2 + \lambda^2 \|\|x\|\|^2 \f$.
	*/
	LSMR& setDamping(const RealScalar& lambda) {
	m_lambda = lambda;
	return *this;
	}

	/** \returns the damping parameter. \sa setDamping() */
	RealScalar damping() const { return m_lambda; }

	/** Sets an upper limit on the estimated condition number of \a A. The
	* iterations stop as soon as the estimate exceeds \a conlim. The default is
	* 0, which disables the limit (equivalent to 1/epsilon).
	*/
	LSMR& setConditionLimit(const RealScalar& conlim) {
	m_conditionLimit = conlim;
	return *this;
	}

	/** \returns the condition-number limit. \sa setConditionLimit() */
	RealScalar conditionLimit() const { return m_conditionLimit; }

	/** Sets the stopping tolerance \c atol, which bounds the relative error
	* assumed in the entries of \a A. It drives the least-squares stopping rule
	* \f$ \|\|A^T r\|\| \le atol\,\|\|A\|\|\,\|\|r\|\| \f$. If left unset (the default) it
	* falls back to tolerance(). \sa setToleranceB(), setTolerance() */
	LSMR& setToleranceA(const RealScalar& atol) {
	m_atol = atol;
	return *this;
	}

	/** \returns \c atol, or tolerance() if setToleranceA() has not been called.
	* \sa setToleranceA() */
	RealScalar toleranceA() const { return m_atol >= RealScalar(0) ? m_atol : Base::m_tolerance; }

	/** Sets the stopping tolerance \c btol, which bounds the relative error
	* assumed in the entries of \a b. It enters the compatible-system stopping
	* rule \f$ \|\|r\|\| \le btol\,\|\|b\|\| + atol\,\|\|A\|\|\,\|\|x\|\| \f$. If left unset (the
	* default) it falls back to tolerance(). \sa setToleranceA(), setTolerance() */
	LSMR& setToleranceB(const RealScalar& btol) {
	m_btol = btol;
	return *this;
	}

	/** \returns \c btol, or tolerance() if setToleranceB() has not been called.
	* \sa setToleranceB() */
	RealScalar toleranceB() const { return m_btol >= RealScalar(0) ? m_btol : Base::m_tolerance; }

	/** \internal */
	template <typename Rhs, typename Dest>
	void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const {
	m_iterations = Base::maxIterations();

	Index istop = internal::lsmr(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error, toleranceA(),
	toleranceB(), m_lambda, m_conditionLimit);
	// istop in {0,1,2,4,5}: the (least-squares) solution was found, possibly
	// only to within machine precision (4,5). istop in {3,6,7}: stopped on the
	// condition-number limit or the iteration limit without meeting the
	// requested tolerance.
	m_info = (istop == 3 \|\| istop == 6 \|\| istop == 7) ? NoConvergence : Success;
	}

	protected:
	RealScalar m_lambda = RealScalar(0);
	RealScalar m_conditionLimit = RealScalar(0);
	// Negative means "unset": toleranceA()/toleranceB() then fall back to tolerance().
	RealScalar m_atol = RealScalar(-1);
	RealScalar m_btol = RealScalar(-1);
	};

	} // end namespace Eigen

	#endif // EIGEN_LSMR_H