Add LSMR iterative least-squares solver

libeigen/eigen!2663

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

Eigen/IterativeLinearSolvers

@@ -27,6 +27,7 @@
   * Those solvers are accessible via the following classes:
   *  - ConjugateGradient for selfadjoint (hermitian) matrices,
   *  - LeastSquaresConjugateGradient for rectangular least-square problems,
+  *  - LSMR for rectangular least-square problems (Golub-Kahan bidiagonalization, optionally damped),
   *  - BiCGSTAB for general square matrices,
   *  - GMRES - a Householder GMRES implementation,
   *  - DGMRES - a deflated GMRES implementation,
@@ -103,6 +104,7 @@
 #include "src/IterativeLinearSolvers/BasicPreconditioners.h"
 #include "src/IterativeLinearSolvers/ConjugateGradient.h"
 #include "src/IterativeLinearSolvers/LeastSquareConjugateGradient.h"
+#include "src/IterativeLinearSolvers/LSMR.h"
 #include "src/IterativeLinearSolvers/BiCGSTAB.h"
 #include "src/IterativeLinearSolvers/IncompleteLUT.h"
 #include "src/IterativeLinearSolvers/IncompleteCholesky.h"

Eigen/src/IterativeLinearSolvers/LSMR.h

@@ -0,0 +1,458 @@
+// 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;
+  typedef typename Dest::RealScalar RealScalar;
+  typedef typename Dest::Scalar Scalar;
+  typedef Matrix<Scalar, Dynamic, 1> VectorType;
+
+  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_> > {
+  typedef MatrixType_ MatrixType;
+  typedef 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().
+ *
+ * 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:
+  typedef IterativeSolverBase<LSMR> Base;
+  using Base::m_error;
+  using Base::m_info;
+  using Base::m_isInitialized;
+  using Base::m_iterations;
+  using Base::matrix;
+
+ public:
+  typedef MatrixType_ MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef 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

doc/SparseLinearSystems.dox

@@ -44,6 +44,10 @@ They are summarized in the following tables:
     <td>IdentityPreconditioner, [LeastSquareDiagonalPreconditioner]</td>
     <td>Solve for min |Ax-b|^2 without forming A'A</td></tr>
 
+<tr><td>LSMR \n <tt>\#include<Eigen/\link IterativeLinearSolvers_Module IterativeLinearSolvers\endlink></tt></td><td>Golub-Kahan bidiagonalization (MINRES on the normal equations)</td><td>Rectangular</td>
+    <td>[IdentityPreconditioner], (self-adjoint right preconditioner)</td>
+    <td>Solve for min |Ax-b|^2, monotonic |A'r|, optional damping. Robust on ill-conditioned problems</td></tr>
+
 <tr><td>BiCGSTAB \n <tt>\#include<Eigen/\link IterativeLinearSolvers_Module IterativeLinearSolvers\endlink></tt></td><td>Iterative stabilized bi-conjugate gradient</td><td>Square</td>
     <td>IdentityPreconditioner, [DiagonalPreconditioner], IncompleteLUT</td>
     <td>To speedup the convergence, try it with the \ref IncompleteLUT preconditioner.</td></tr>

setup.cfg

@@ -15,6 +15,7 @@ skip = .git,build,.build,build-clang,*.o,*.a,*.so,*.pyc,blas/f2c,blas/testing,la
 # Tensor parameter name: vaccum (vector accumulator)
 # Data-structure plural: deques (correct plural of deque)
 # Cholesky method name: applyD (apply the block-diagonal matrix D)
+# LSMR damping-rotation sine: shat (s-hat, from Fong & Saunders)
 ignore-words-list =
     nd,te,ba,ser,
     mata,matc,
@@ -28,4 +29,5 @@ ignore-words-list =
     bload,vaccum,
     deques,
     applyd,
+    shat,
     manuel

test/CMakeLists.txt

@@ -324,6 +324,7 @@ ei_add_test(incomplete_cholesky)
 ei_add_test(incomplete_LUT)
 ei_add_test(bicgstab)
 ei_add_test(lscg)
+ei_add_test(lsmr)
 ei_add_test(gmres)
 ei_add_test(dgmres)
 ei_add_test(minres)

test/lsmr.cpp

@@ -0,0 +1,240 @@
+// 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
+
+#include "sparse_solver.h"
+#include <Eigen/IterativeLinearSolvers>
+
+// Sparse square and least-squares solving with the default (identity) and the
+// least-squares diagonal preconditioner, in both storage orders.
+template <typename T>
+void test_lsmr_T() {
+  LSMR<SparseMatrix<T> > lsmr_colmajor_I;
+  LSMR<SparseMatrix<T>, LeastSquareDiagonalPreconditioner<T> > lsmr_colmajor_diag;
+  LSMR<SparseMatrix<T, RowMajor> > lsmr_rowmajor_I;
+  LSMR<SparseMatrix<T, RowMajor>, LeastSquareDiagonalPreconditioner<T> > lsmr_rowmajor_diag;
+
+  CALL_SUBTEST(check_sparse_square_solving(lsmr_colmajor_I));
+  CALL_SUBTEST(check_sparse_square_solving(lsmr_colmajor_diag));
+
+  CALL_SUBTEST(check_sparse_leastsquare_solving(lsmr_colmajor_I));
+  CALL_SUBTEST(check_sparse_leastsquare_solving(lsmr_colmajor_diag));
+
+  CALL_SUBTEST(check_sparse_square_solving(lsmr_rowmajor_I));
+  CALL_SUBTEST(check_sparse_leastsquare_solving(lsmr_rowmajor_diag));
+}
+
+// LSMR on dense, rectangular, overdetermined problems: the solution must match
+// the dense QR least-squares solution.
+template <typename T>
+void test_lsmr_dense() {
+  typedef typename NumTraits<T>::Real RealScalar;
+  typedef Matrix<T, Dynamic, Dynamic> DenseMatrix;
+  typedef Matrix<T, Dynamic, 1> DenseVector;
+
+  for (int k = 0; k < g_repeat; ++k) {
+    Index rows = internal::random<Index>(20, 80);
+    Index cols = internal::random<Index>(4, rows);
+    DenseMatrix A = DenseMatrix::Random(rows, cols);
+    DenseVector b = DenseVector::Random(rows);
+
+    LSMR<DenseMatrix> lsmr(A);
+    lsmr.setTolerance(NumTraits<RealScalar>::epsilon() * RealScalar(100));
+    DenseVector x = lsmr.solve(b);
+    VERIFY_IS_EQUAL(lsmr.info(), Success);
+
+    DenseVector xref = A.householderQr().solve(b);
+    VERIFY_IS_APPROX(x, xref);
+  }
+}
+
+// solveWithGuess() must actually consume the initial guess. Starting from the
+// exact solution, the normal-equation residual A^T(b - A x0) is ~0, so LSMR
+// converges in very few iterations -- far fewer than from x = 0. Were the guess
+// silently dropped, both runs would start from x = 0 and use the same count, so
+// the strict iteration-count comparison below is load-bearing. We use cols >= 10
+// so that the from-zero run needs many iterations, giving a wide margin.
+template <typename T>
+void test_lsmr_guess() {
+  typedef typename NumTraits<T>::Real RealScalar;
+  typedef Matrix<T, Dynamic, Dynamic> DenseMatrix;
+  typedef Matrix<T, Dynamic, 1> DenseVector;
+
+  for (int k = 0; k < g_repeat; ++k) {
+    Index cols = internal::random<Index>(10, 20);
+    Index rows = cols + internal::random<Index>(10, 40);
+    DenseMatrix A = DenseMatrix::Random(rows, cols);
+    DenseVector b = DenseVector::Random(rows);
+    DenseVector xref = A.householderQr().solve(b);
+
+    LSMR<DenseMatrix> lsmr(A);
+    lsmr.setTolerance(NumTraits<RealScalar>::epsilon() * RealScalar(100));
+
+    DenseVector x0 = lsmr.solve(b);
+    VERIFY_IS_EQUAL(lsmr.info(), Success);
+    VERIFY_IS_APPROX(x0, xref);
+    const Index iters_from_zero = lsmr.iterations();
+
+    DenseVector xg = lsmr.solveWithGuess(b, xref);
+    VERIFY_IS_EQUAL(lsmr.info(), Success);
+    VERIFY_IS_APPROX(xg, xref);
+    VERIFY(lsmr.iterations() < iters_from_zero);
+  }
+}
+
+// On rank-deficient or underdetermined (rows < cols) problems the least-squares
+// solution is not unique. Started from x = 0, LSMR returns the *minimum-norm*
+// least-squares solution -- exactly what completeOrthogonalDecomposition
+// returns, whereas householderQr would return some other valid LS solution.
+template <typename T>
+void test_lsmr_minnorm() {
+  typedef typename NumTraits<T>::Real RealScalar;
+  typedef Matrix<T, Dynamic, Dynamic> DenseMatrix;
+  typedef Matrix<T, Dynamic, 1> DenseVector;
+
+  for (int k = 0; k < g_repeat; ++k) {
+    // (a) Underdetermined, full row rank: rows < cols, many exact solutions.
+    {
+      Index rows = internal::random<Index>(10, 30);
+      Index cols = rows + internal::random<Index>(5, 30);
+      DenseMatrix A = DenseMatrix::Random(rows, cols);
+      DenseVector b = DenseVector::Random(rows);
+
+      LSMR<DenseMatrix> lsmr(A);
+      lsmr.setTolerance(NumTraits<RealScalar>::epsilon() * RealScalar(100));
+      DenseVector x = lsmr.solve(b);
+      VERIFY_IS_EQUAL(lsmr.info(), Success);
+
+      DenseVector xref = A.completeOrthogonalDecomposition().solve(b);
+      VERIFY_IS_APPROX(x, xref);
+    }
+    // (b) Explicitly rank-deficient, overdetermined: a duplicated column.
+    {
+      Index rows = internal::random<Index>(30, 60);
+      Index cols = internal::random<Index>(5, 15);
+      DenseMatrix A = DenseMatrix::Random(rows, cols);
+      A.col(1) = A.col(0);  // make A rank-deficient
+      DenseVector b = DenseVector::Random(rows);
+
+      LSMR<DenseMatrix> lsmr(A);
+      lsmr.setTolerance(NumTraits<RealScalar>::epsilon() * RealScalar(100));
+      DenseVector x = lsmr.solve(b);
+      VERIFY_IS_EQUAL(lsmr.info(), Success);
+
+      DenseVector xref = A.completeOrthogonalDecomposition().solve(b);
+      VERIFY_IS_APPROX(x, xref);
+    }
+  }
+}
+
+// With damping lambda > 0 LSMR solves the regularized problem
+//   min ||Ax-b||^2 + lambda^2 ||x||^2,
+// whose solution satisfies the (well-conditioned) normal equation
+//   (A^T A + lambda^2 I) x = A^T b.
+template <typename T>
+void test_lsmr_damping() {
+  typedef typename NumTraits<T>::Real RealScalar;
+  typedef Matrix<T, Dynamic, Dynamic> DenseMatrix;
+  typedef Matrix<T, Dynamic, 1> DenseVector;
+
+  for (int k = 0; k < g_repeat; ++k) {
+    Index rows = internal::random<Index>(20, 60);
+    Index cols = internal::random<Index>(4, rows);
+    DenseMatrix A = DenseMatrix::Random(rows, cols);
+    DenseVector b = DenseVector::Random(rows);
+    RealScalar lambda = internal::random<RealScalar>(RealScalar(0.2), RealScalar(2));
+
+    LSMR<DenseMatrix> lsmr(A);
+    lsmr.setDamping(lambda);
+    lsmr.setTolerance(NumTraits<RealScalar>::epsilon() * RealScalar(100));
+    DenseVector x = lsmr.solve(b);
+    VERIFY_IS_EQUAL(lsmr.info(), Success);
+
+    DenseMatrix normalMat = A.adjoint() * A + (lambda * lambda) * DenseMatrix::Identity(cols, cols);
+    DenseVector xref = normalMat.ldlt().solve(A.adjoint() * b);
+    VERIFY_IS_APPROX(x, xref);
+  }
+}
+
+// Robustness against extreme right-hand-side scaling: LSMR normalizes by ||b||
+// immediately, so the solution should scale linearly with b.
+void test_lsmr_extreme_rhs() {
+  const Matrix2d mat = Matrix2d::Identity();
+  const Vector2d direction = (Vector2d() << 1, -1).finished();
+  LSMR<Matrix2d> solver(mat);
+  solver.setTolerance(1e-12);
+
+  for (double scale : {1e-200, 1e200}) {
+    const Vector2d rhs = scale * direction;
+    Vector2d x = solver.solve(rhs);
+    VERIFY_IS_EQUAL(solver.info(), Success);
+    VERIFY(x.allFinite());
+    VERIFY_IS_APPROX(x / scale, direction);
+  }
+}
+
+// The atol/btol stopping tolerances can be set independently. setTolerance()
+// sets both; setToleranceA()/setToleranceB() override each one, and an unset
+// component falls back to tolerance(). A tight atol must still reach the QR
+// least-squares solution even when btol is left loose.
+template <typename T>
+void test_lsmr_tolerances() {
+  typedef typename NumTraits<T>::Real RealScalar;
+  typedef Matrix<T, Dynamic, Dynamic> DenseMatrix;
+  typedef Matrix<T, Dynamic, 1> DenseVector;
+
+  const RealScalar def = NumTraits<RealScalar>::epsilon();
+  const RealScalar tight = def * RealScalar(100);
+
+  // By default both tolerances track tolerance().
+  LSMR<DenseMatrix> lsmr;
+  VERIFY_IS_EQUAL(lsmr.toleranceA(), def);
+  VERIFY_IS_EQUAL(lsmr.toleranceB(), def);
+
+  // setTolerance() sets both.
+  lsmr.setTolerance(tight);
+  VERIFY_IS_EQUAL(lsmr.toleranceA(), tight);
+  VERIFY_IS_EQUAL(lsmr.toleranceB(), tight);
+
+  // An explicit override wins; the other component still tracks tolerance().
+  lsmr.setToleranceA(def);
+  VERIFY_IS_EQUAL(lsmr.toleranceA(), def);
+  VERIFY_IS_EQUAL(lsmr.toleranceB(), tight);
+  lsmr.setToleranceB(def);
+  VERIFY_IS_EQUAL(lsmr.toleranceB(), def);
+
+  for (int k = 0; k < g_repeat; ++k) {
+    Index rows = internal::random<Index>(20, 80);
+    Index cols = internal::random<Index>(4, rows);
+    DenseMatrix A = DenseMatrix::Random(rows, cols);
+    DenseVector b = DenseVector::Random(rows);
+    DenseVector xref = A.householderQr().solve(b);
+
+    // Tight atol drives the least-squares stopping rule, so the solution must
+    // match QR even with btol left loose.
+    LSMR<DenseMatrix> solver(A);
+    solver.setToleranceA(tight).setToleranceB(RealScalar(0.1));
+    DenseVector x = solver.solve(b);
+    VERIFY_IS_EQUAL(solver.info(), Success);
+    VERIFY_IS_APPROX(x, xref);
+  }
+}
+
+EIGEN_DECLARE_TEST(lsmr) {
+  CALL_SUBTEST_1(test_lsmr_T<double>());
+  CALL_SUBTEST_2(test_lsmr_T<std::complex<double> >());
+  CALL_SUBTEST_3(test_lsmr_dense<double>());
+  CALL_SUBTEST_4(test_lsmr_dense<std::complex<double> >());
+  CALL_SUBTEST_5(test_lsmr_damping<double>());
+  CALL_SUBTEST_6(test_lsmr_extreme_rhs());
+  CALL_SUBTEST_7(test_lsmr_minnorm<double>());
+  CALL_SUBTEST_8(test_lsmr_damping<std::complex<double> >());
+  CALL_SUBTEST_9(test_lsmr_minnorm<std::complex<double> >());
+  CALL_SUBTEST_10(test_lsmr_guess<double>());
+  CALL_SUBTEST_11(test_lsmr_tolerances<double>());
+}