src/IterativeSolvers/MINRES.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Giacomo Po <gpo@ucla.edu>
0005 // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 // Copyright (C) 2018 David Hyde <dabh@stanford.edu>
0007 //
0008 // This Source Code Form is subject to the terms of the Mozilla
0009 // Public License v. 2.0. If a copy of the MPL was not distributed
0010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0011
0012
0013 #ifndef EIGEN_MINRES_H_
0014 #define EIGEN_MINRES_H_
0015
0016
0017 namespace Eigen {
0018
0019     namespace internal {
0020
0021         /** \internal Low-level MINRES algorithm
0022          * \param mat The matrix A
0023          * \param rhs The right hand side vector b
0024          * \param x On input and initial solution, on output the computed solution.
0025          * \param precond A right preconditioner being able to efficiently solve for an
0026          *                approximation of Ax=b (regardless of b)
0027          * \param iters On input the max number of iteration, on output the number of performed iterations.
0028          * \param tol_error On input the tolerance error, on output an estimation of the relative error.
0029          */
0030         template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
0031         EIGEN_DONT_INLINE
0032         void minres(const MatrixType& mat, const Rhs& rhs, Dest& x,
0033                     const Preconditioner& precond, Index& iters,
0034                     typename Dest::RealScalar& tol_error)
0035         {
0036             using std::sqrt;
0037             typedef typename Dest::RealScalar RealScalar;
0038             typedef typename Dest::Scalar Scalar;
0039             typedef Matrix<Scalar,Dynamic,1> VectorType;
0040
0041             // Check for zero rhs
0042             const RealScalar rhsNorm2(rhs.squaredNorm());
0043             if(rhsNorm2 == 0)
0044             {
0045                 x.setZero();
0046                 iters = 0;
0047                 tol_error = 0;
0048                 return;
0049             }
0050
0051             // initialize
0052             const Index maxIters(iters);  // initialize maxIters to iters
0053             const Index N(mat.cols());    // the size of the matrix
0054             const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2)
0055
0056             // Initialize preconditioned Lanczos
0057             VectorType v_old(N); // will be initialized inside loop
0058             VectorType v( VectorType::Zero(N) ); //initialize v
0059             VectorType v_new(rhs-mat*x); //initialize v_new
0060             RealScalar residualNorm2(v_new.squaredNorm());
0061             VectorType w(N); // will be initialized inside loop
0062             VectorType w_new(precond.solve(v_new)); // initialize w_new
0063 //            RealScalar beta; // will be initialized inside loop
0064             RealScalar beta_new2(v_new.dot(w_new));
0065             eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
0066             RealScalar beta_new(sqrt(beta_new2));
0067             const RealScalar beta_one(beta_new);
0068             // Initialize other variables
0069             RealScalar c(1.0); // the cosine of the Givens rotation
0070             RealScalar c_old(1.0);
0071             RealScalar s(0.0); // the sine of the Givens rotation
0072             RealScalar s_old(0.0); // the sine of the Givens rotation
0073             VectorType p_oold(N); // will be initialized in loop
0074             VectorType p_old(VectorType::Zero(N)); // initialize p_old=0
0075             VectorType p(p_old); // initialize p=0
0076             RealScalar eta(1.0);
0077
0078             iters = 0; // reset iters
0079             while ( iters < maxIters )
0080             {
0081                 // Preconditioned Lanczos
0082                 /* Note that there are 4 variants on the Lanczos algorithm. These are
0083                  * described in Paige, C. C. (1972). Computational variants of
0084                  * the Lanczos method for the eigenproblem. IMA Journal of Applied
0085                  * Mathematics, 10(3), 373-381. The current implementation corresponds
0086                  * to the case A(2,7) in the paper. It also corresponds to
0087                  * algorithm 6.14 in Y. Saad, Iterative Methods for Sparse Linear
0088                  * Systems, 2003 p.173. For the preconditioned version see
0089                  * A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987).
0090                  */
0091                 const RealScalar beta(beta_new);
0092                 v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter
0093                 v_new /= beta_new; // overwrite v_new for next iteration
0094                 w_new /= beta_new; // overwrite w_new for next iteration
0095                 v = v_new; // update
0096                 w = w_new; // update
0097                 v_new.noalias() = mat*w - beta*v_old; // compute v_new
0098                 const RealScalar alpha = v_new.dot(w);
0099                 v_new -= alpha*v; // overwrite v_new
0100                 w_new = precond.solve(v_new); // overwrite w_new
0101                 beta_new2 = v_new.dot(w_new); // compute beta_new
0102                 eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
0103                 beta_new = sqrt(beta_new2); // compute beta_new
0104
0105                 // Givens rotation
0106                 const RealScalar r2 =s*alpha+c*c_old*beta; // s, s_old, c and c_old are still from previous iteration
0107                 const RealScalar r3 =s_old*beta; // s, s_old, c and c_old are still from previous iteration
0108                 const RealScalar r1_hat=c*alpha-c_old*s*beta;
0109                 const RealScalar r1 =sqrt( std::pow(r1_hat,2) + std::pow(beta_new,2) );
0110                 c_old = c; // store for next iteration
0111                 s_old = s; // store for next iteration
0112                 c=r1_hat/r1; // new cosine
0113                 s=beta_new/r1; // new sine
0114
0115                 // Update solution
0116                 p_oold = p_old;
0117                 p_old = p;
0118                 p.noalias()=(w-r2*p_old-r3*p_oold) /r1; // IS NOALIAS REQUIRED?
0119                 x += beta_one*c*eta*p;
0120
0121                 /* Update the squared residual. Note that this is the estimated residual.
0122                 The real residual |Ax-b|^2 may be slightly larger */
0123                 residualNorm2 *= s*s;
0124
0125                 if ( residualNorm2 < threshold2)
0126                 {
0127                     break;
0128                 }
0129
0130                 eta=-s*eta; // update eta
0131                 iters++; // increment iteration number (for output purposes)
0132             }
0133
0134             /* Compute error. Note that this is the estimated error. The real
0135              error |Ax-b|/|b| may be slightly larger */
0136             tol_error = std::sqrt(residualNorm2 / rhsNorm2);
0137         }
0138
0139     }
0140
0141     template< typename _MatrixType, int _UpLo=Lower,
0142     typename _Preconditioner = IdentityPreconditioner>
0143     class MINRES;
0144
0145     namespace internal {
0146
0147         template< typename _MatrixType, int _UpLo, typename _Preconditioner>
0148         struct traits<MINRES<_MatrixType,_UpLo,_Preconditioner> >
0149         {
0150             typedef _MatrixType MatrixType;
0151             typedef _Preconditioner Preconditioner;
0152         };
0153
0154     }
0155
0156     /** \ingroup IterativeLinearSolvers_Module
0157      * \brief A minimal residual solver for sparse symmetric problems
0158      *
0159      * This class allows to solve for A.x = b sparse linear problems using the MINRES algorithm
0160      * of Paige and Saunders (1975). The sparse matrix A must be symmetric (possibly indefinite).
0161      * The vectors x and b can be either dense or sparse.
0162      *
0163      * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix.
0164      * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower,
0165      *               Upper, or Lower|Upper in which the full matrix entries will be considered. Default is Lower.
0166      * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
0167      *
0168      * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
0169      * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
0170      * and NumTraits<Scalar>::epsilon() for the tolerance.
0171      *
0172      * This class can be used as the direct solver classes. Here is a typical usage example:
0173      * \code
0174      * int n = 10000;
0175      * VectorXd x(n), b(n);
0176      * SparseMatrix<double> A(n,n);
0177      * // fill A and b
0178      * MINRES<SparseMatrix<double> > mr;
0179      * mr.compute(A);
0180      * x = mr.solve(b);
0181      * std::cout << "#iterations:     " << mr.iterations() << std::endl;
0182      * std::cout << "estimated error: " << mr.error()      << std::endl;
0183      * // update b, and solve again
0184      * x = mr.solve(b);
0185      * \endcode
0186      *
0187      * By default the iterations start with x=0 as an initial guess of the solution.
0188      * One can control the start using the solveWithGuess() method.
0189      *
0190      * MINRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
0191      *
0192      * \sa class ConjugateGradient, BiCGSTAB, SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
0193      */
0194     template< typename _MatrixType, int _UpLo, typename _Preconditioner>
0195     class MINRES : public IterativeSolverBase<MINRES<_MatrixType,_UpLo,_Preconditioner> >
0196     {
0197
0198         typedef IterativeSolverBase<MINRES> Base;
0199         using Base::matrix;
0200         using Base::m_error;
0201         using Base::m_iterations;
0202         using Base::m_info;
0203         using Base::m_isInitialized;
0204     public:
0205         using Base::_solve_impl;
0206         typedef _MatrixType MatrixType;
0207         typedef typename MatrixType::Scalar Scalar;
0208         typedef typename MatrixType::RealScalar RealScalar;
0209         typedef _Preconditioner Preconditioner;
0210
0211         enum {UpLo = _UpLo};
0212
0213     public:
0214
0215         /** Default constructor. */
0216         MINRES() : Base() {}
0217
0218         /** Initialize the solver with matrix \a A for further \c Ax=b solving.
0219          *
0220          * This constructor is a shortcut for the default constructor followed
0221          * by a call to compute().
0222          *
0223          * \warning this class stores a reference to the matrix A as well as some
0224          * precomputed values that depend on it. Therefore, if \a A is changed
0225          * this class becomes invalid. Call compute() to update it with the new
0226          * matrix A, or modify a copy of A.
0227          */
0228         template<typename MatrixDerived>
0229         explicit MINRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
0230
0231         /** Destructor. */
0232         ~MINRES(){}
0233
0234         /** \internal */
0235         template<typename Rhs,typename Dest>
0236         void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
0237         {
0238             typedef typename Base::MatrixWrapper MatrixWrapper;
0239             typedef typename Base::ActualMatrixType ActualMatrixType;
0240             enum {
0241               TransposeInput  =   (!MatrixWrapper::MatrixFree)
0242                               &&  (UpLo==(Lower|Upper))
0243                               &&  (!MatrixType::IsRowMajor)
0244                               &&  (!NumTraits<Scalar>::IsComplex)
0245             };
0246             typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType const&>::type RowMajorWrapper;
0247             EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
0248             typedef typename internal::conditional<UpLo==(Lower|Upper),
0249                                                   RowMajorWrapper,
0250                                                   typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
0251                                             >::type SelfAdjointWrapper;
0252
0253             m_iterations = Base::maxIterations();
0254             m_error = Base::m_tolerance;
0255             RowMajorWrapper row_mat(matrix());
0256             internal::minres(SelfAdjointWrapper(row_mat), b, x,
0257                              Base::m_preconditioner, m_iterations, m_error);
0258             m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
0259         }
0260
0261     protected:
0262
0263     };
0264
0265 } // end namespace Eigen
0266
0267 #endif // EIGEN_MINRES_H