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 {

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

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_;

  LSMR() : Base() {}

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

  LSMR& setDamping(const RealScalar& lambda) {
    m_lambda = lambda;
    return *this;
  }

  RealScalar damping() const { return m_lambda; }

  LSMR& setConditionLimit(const RealScalar& conlim) {
    m_conditionLimit = conlim;
    return *this;
  }

  RealScalar conditionLimit() const { return m_conditionLimit; }

  LSMR& setToleranceA(const RealScalar& atol) {
    m_atol = atol;
    return *this;
  }

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

  LSMR& setToleranceB(const RealScalar& btol) {
    m_btol = btol;
    return *this;
  }

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

  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