1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2015 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_LEAST_SQUARE_CONJUGATE_GRADIENT_H 11 #define EIGEN_LEAST_SQUARE_CONJUGATE_GRADIENT_H 12 13 namespace Eigen { 14 15 namespace internal { 16 17 /** \internal Low-level conjugate gradient algorithm for least-square problems 18 * \param mat The matrix A 19 * \param rhs The right hand side vector b 20 * \param x On input and initial solution, on output the computed solution. 21 * \param precond A preconditioner being able to efficiently solve for an 22 * approximation of A'Ax=b (regardless of b) 23 * \param iters On input the max number of iteration, on output the number of performed iterations. 24 * \param tol_error On input the tolerance error, on output an estimation of the relative error. 25 */ 26 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> 27 EIGEN_DONT_INLINE 28 void least_square_conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, 29 const Preconditioner& precond, Index& iters, 30 typename Dest::RealScalar& tol_error) 31 { 32 using std::sqrt; 33 using std::abs; 34 typedef typename Dest::RealScalar RealScalar; 35 typedef typename Dest::Scalar Scalar; 36 typedef Matrix<Scalar,Dynamic,1> VectorType; 37 38 RealScalar tol = tol_error; 39 Index maxIters = iters; 40 41 Index m = mat.rows(), n = mat.cols(); 42 43 VectorType residual = rhs - mat * x; 44 VectorType normal_residual = mat.adjoint() * residual; 45 46 RealScalar rhsNorm2 = (mat.adjoint()*rhs).squaredNorm(); 47 if(rhsNorm2 == 0) 48 { 49 x.setZero(); 50 iters = 0; 51 tol_error = 0; 52 return; 53 } 54 RealScalar threshold = tol*tol*rhsNorm2; 55 RealScalar residualNorm2 = normal_residual.squaredNorm(); 56 if (residualNorm2 < threshold) 57 { 58 iters = 0; 59 tol_error = sqrt(residualNorm2 / rhsNorm2); 60 return; 61 } 62 63 VectorType p(n); 64 p = precond.solve(normal_residual); // initial search direction 65 66 VectorType z(n), tmp(m); 67 RealScalar absNew = numext::real(normal_residual.dot(p)); // the square of the absolute value of r scaled by invM 68 Index i = 0; 69 while(i < maxIters) 70 { 71 tmp.noalias() = mat * p; 72 73 Scalar alpha = absNew / tmp.squaredNorm(); // the amount we travel on dir 74 x += alpha * p; // update solution 75 residual -= alpha * tmp; // update residual 76 normal_residual = mat.adjoint() * residual; // update residual of the normal equation 77 78 residualNorm2 = normal_residual.squaredNorm(); 79 if(residualNorm2 < threshold) 80 break; 81 82 z = precond.solve(normal_residual); // approximately solve for "A'A z = normal_residual" 83 84 RealScalar absOld = absNew; 85 absNew = numext::real(normal_residual.dot(z)); // update the absolute value of r 86 RealScalar beta = absNew / absOld; // calculate the Gram-Schmidt value used to create the new search direction 87 p = z + beta * p; // update search direction 88 i++; 89 } 90 tol_error = sqrt(residualNorm2 / rhsNorm2); 91 iters = i; 92 } 93 94 } 95 96 template< typename _MatrixType, 97 typename _Preconditioner = LeastSquareDiagonalPreconditioner<typename _MatrixType::Scalar> > 98 class LeastSquaresConjugateGradient; 99 100 namespace internal { 101 102 template< typename _MatrixType, typename _Preconditioner> 103 struct traits<LeastSquaresConjugateGradient<_MatrixType,_Preconditioner> > 104 { 105 typedef _MatrixType MatrixType; 106 typedef _Preconditioner Preconditioner; 107 }; 108 109 } 110 111 /** \ingroup IterativeLinearSolvers_Module 112 * \brief A conjugate gradient solver for sparse (or dense) least-square problems 113 * 114 * This class allows to solve for A x = b linear problems using an iterative conjugate gradient algorithm. 115 * The matrix A can be non symmetric and rectangular, but the matrix A' A should be positive-definite to guaranty stability. 116 * Otherwise, the SparseLU or SparseQR classes might be preferable. 117 * The matrix A and the vectors x and b can be either dense or sparse. 118 * 119 * \tparam _MatrixType the type of the matrix A, can be a dense or a sparse matrix. 120 * \tparam _Preconditioner the type of the preconditioner. Default is LeastSquareDiagonalPreconditioner 121 * 122 * \implsparsesolverconcept 123 * 124 * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() 125 * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations 126 * and NumTraits<Scalar>::epsilon() for the tolerance. 127 * 128 * This class can be used as the direct solver classes. Here is a typical usage example: 129 \code 130 int m=1000000, n = 10000; 131 VectorXd x(n), b(m); 132 SparseMatrix<double> A(m,n); 133 // fill A and b 134 LeastSquaresConjugateGradient<SparseMatrix<double> > lscg; 135 lscg.compute(A); 136 x = lscg.solve(b); 137 std::cout << "#iterations: " << lscg.iterations() << std::endl; 138 std::cout << "estimated error: " << lscg.error() << std::endl; 139 // update b, and solve again 140 x = lscg.solve(b); 141 \endcode 142 * 143 * By default the iterations start with x=0 as an initial guess of the solution. 144 * One can control the start using the solveWithGuess() method. 145 * 146 * \sa class ConjugateGradient, SparseLU, SparseQR 147 */ 148 template< typename _MatrixType, typename _Preconditioner> 149 class LeastSquaresConjugateGradient : public IterativeSolverBase<LeastSquaresConjugateGradient<_MatrixType,_Preconditioner> > 150 { 151 typedef IterativeSolverBase<LeastSquaresConjugateGradient> Base; 152 using Base::matrix; 153 using Base::m_error; 154 using Base::m_iterations; 155 using Base::m_info; 156 using Base::m_isInitialized; 157 public: 158 typedef _MatrixType MatrixType; 159 typedef typename MatrixType::Scalar Scalar; 160 typedef typename MatrixType::RealScalar RealScalar; 161 typedef _Preconditioner Preconditioner; 162 163 public: 164 165 /** Default constructor. */ 166 LeastSquaresConjugateGradient() : Base() {} 167 168 /** Initialize the solver with matrix \a A for further \c Ax=b solving. 169 * 170 * This constructor is a shortcut for the default constructor followed 171 * by a call to compute(). 172 * 173 * \warning this class stores a reference to the matrix A as well as some 174 * precomputed values that depend on it. Therefore, if \a A is changed 175 * this class becomes invalid. Call compute() to update it with the new 176 * matrix A, or modify a copy of A. 177 */ 178 template<typename MatrixDerived> 179 explicit LeastSquaresConjugateGradient(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {} 180 181 ~LeastSquaresConjugateGradient() {} 182 183 /** \internal */ 184 template<typename Rhs,typename Dest> 185 void _solve_with_guess_impl(const Rhs& b, Dest& x) const 186 { 187 m_iterations = Base::maxIterations(); 188 m_error = Base::m_tolerance; 189 190 for(Index j=0; j<b.cols(); ++j) 191 { 192 m_iterations = Base::maxIterations(); 193 m_error = Base::m_tolerance; 194 195 typename Dest::ColXpr xj(x,j); 196 internal::least_square_conjugate_gradient(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error); 197 } 198 199 m_isInitialized = true; 200 m_info = m_error <= Base::m_tolerance ? Success : NoConvergence; 201 } 202 203 /** \internal */ 204 using Base::_solve_impl; 205 template<typename Rhs,typename Dest> 206 void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const 207 { 208 x.setZero(); 209 _solve_with_guess_impl(b.derived(),x); 210 } 211 212 }; 213 214 } // end namespace Eigen 215 216 #endif // EIGEN_LEAST_SQUARE_CONJUGATE_GRADIENT_H 217
