diff --git a/Eigen/IterativeLinearSolvers b/Eigen/IterativeLinearSolvers index 882c49d52..15631ef41 100644 --- a/Eigen/IterativeLinearSolvers +++ b/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" diff --git a/Eigen/src/IterativeLinearSolvers/LSMR.h b/Eigen/src/IterativeLinearSolvers/LSMR.h new file mode 100644 index 000000000..b328130f5 --- /dev/null +++ b/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 +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 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::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 +class LSMR; + +namespace internal { + +template +struct traits > { + 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::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 A(m, n); + // fill A and b + LSMR > 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 +class LSMR : public IterativeSolverBase > { + protected: + typedef IterativeSolverBase 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 + explicit LSMR(const EigenBase& 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 + 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 diff --git a/doc/SparseLinearSystems.dox b/doc/SparseLinearSystems.dox index 0c7bfc15c..ed6e1b5cd 100644 --- a/doc/SparseLinearSystems.dox +++ b/doc/SparseLinearSystems.dox @@ -44,6 +44,10 @@ They are summarized in the following tables: IdentityPreconditioner, [LeastSquareDiagonalPreconditioner] Solve for min |Ax-b|^2 without forming A'A +LSMR \n \#includeGolub-Kahan bidiagonalization (MINRES on the normal equations)Rectangular + [IdentityPreconditioner], (self-adjoint right preconditioner) + Solve for min |Ax-b|^2, monotonic |A'r|, optional damping. Robust on ill-conditioned problems + BiCGSTAB \n \#includeIterative stabilized bi-conjugate gradientSquare IdentityPreconditioner, [DiagonalPreconditioner], IncompleteLUT To speedup the convergence, try it with the \ref IncompleteLUT preconditioner. diff --git a/setup.cfg b/setup.cfg index ced1cc803..960d3b3d0 100644 --- a/setup.cfg +++ b/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 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c4ad4478c..80ed9eb35 100644 --- a/test/CMakeLists.txt +++ b/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) diff --git a/test/lsmr.cpp b/test/lsmr.cpp new file mode 100644 index 000000000..6d9157fd2 --- /dev/null +++ b/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 + +// Sparse square and least-squares solving with the default (identity) and the +// least-squares diagonal preconditioner, in both storage orders. +template +void test_lsmr_T() { + LSMR > lsmr_colmajor_I; + LSMR, LeastSquareDiagonalPreconditioner > lsmr_colmajor_diag; + LSMR > lsmr_rowmajor_I; + LSMR, LeastSquareDiagonalPreconditioner > 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 +void test_lsmr_dense() { + typedef typename NumTraits::Real RealScalar; + typedef Matrix DenseMatrix; + typedef Matrix DenseVector; + + for (int k = 0; k < g_repeat; ++k) { + Index rows = internal::random(20, 80); + Index cols = internal::random(4, rows); + DenseMatrix A = DenseMatrix::Random(rows, cols); + DenseVector b = DenseVector::Random(rows); + + LSMR lsmr(A); + lsmr.setTolerance(NumTraits::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 +void test_lsmr_guess() { + typedef typename NumTraits::Real RealScalar; + typedef Matrix DenseMatrix; + typedef Matrix DenseVector; + + for (int k = 0; k < g_repeat; ++k) { + Index cols = internal::random(10, 20); + Index rows = cols + internal::random(10, 40); + DenseMatrix A = DenseMatrix::Random(rows, cols); + DenseVector b = DenseVector::Random(rows); + DenseVector xref = A.householderQr().solve(b); + + LSMR lsmr(A); + lsmr.setTolerance(NumTraits::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 +void test_lsmr_minnorm() { + typedef typename NumTraits::Real RealScalar; + typedef Matrix DenseMatrix; + typedef Matrix DenseVector; + + for (int k = 0; k < g_repeat; ++k) { + // (a) Underdetermined, full row rank: rows < cols, many exact solutions. + { + Index rows = internal::random(10, 30); + Index cols = rows + internal::random(5, 30); + DenseMatrix A = DenseMatrix::Random(rows, cols); + DenseVector b = DenseVector::Random(rows); + + LSMR lsmr(A); + lsmr.setTolerance(NumTraits::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(30, 60); + Index cols = internal::random(5, 15); + DenseMatrix A = DenseMatrix::Random(rows, cols); + A.col(1) = A.col(0); // make A rank-deficient + DenseVector b = DenseVector::Random(rows); + + LSMR lsmr(A); + lsmr.setTolerance(NumTraits::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 +void test_lsmr_damping() { + typedef typename NumTraits::Real RealScalar; + typedef Matrix DenseMatrix; + typedef Matrix DenseVector; + + for (int k = 0; k < g_repeat; ++k) { + Index rows = internal::random(20, 60); + Index cols = internal::random(4, rows); + DenseMatrix A = DenseMatrix::Random(rows, cols); + DenseVector b = DenseVector::Random(rows); + RealScalar lambda = internal::random(RealScalar(0.2), RealScalar(2)); + + LSMR lsmr(A); + lsmr.setDamping(lambda); + lsmr.setTolerance(NumTraits::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 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 +void test_lsmr_tolerances() { + typedef typename NumTraits::Real RealScalar; + typedef Matrix DenseMatrix; + typedef Matrix DenseVector; + + const RealScalar def = NumTraits::epsilon(); + const RealScalar tight = def * RealScalar(100); + + // By default both tolerances track tolerance(). + LSMR 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(20, 80); + Index cols = internal::random(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 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()); + CALL_SUBTEST_2(test_lsmr_T >()); + CALL_SUBTEST_3(test_lsmr_dense()); + CALL_SUBTEST_4(test_lsmr_dense >()); + CALL_SUBTEST_5(test_lsmr_damping()); + CALL_SUBTEST_6(test_lsmr_extreme_rhs()); + CALL_SUBTEST_7(test_lsmr_minnorm()); + CALL_SUBTEST_8(test_lsmr_damping >()); + CALL_SUBTEST_9(test_lsmr_minnorm >()); + CALL_SUBTEST_10(test_lsmr_guess()); + CALL_SUBTEST_11(test_lsmr_tolerances()); +}