Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
LeastSquareConjugateGradient.h
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@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
13namespace Eigen {
14
15namespace internal {
16
26template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
27EIGEN_DONT_INLINE
28void 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
96template< typename _MatrixType,
97 typename _Preconditioner = LeastSquareDiagonalPreconditioner<typename _MatrixType::Scalar> >
98class LeastSquaresConjugateGradient;
99
100namespace internal {
101
102template< typename _MatrixType, typename _Preconditioner>
103struct traits<LeastSquaresConjugateGradient<_MatrixType,_Preconditioner> >
104{
105 typedef _MatrixType MatrixType;
106 typedef _Preconditioner Preconditioner;
107};
108
109}
110
148template< typename _MatrixType, typename _Preconditioner>
149class LeastSquaresConjugateGradient : public IterativeSolverBase<LeastSquaresConjugateGradient<_MatrixType,_Preconditioner> >
150{
152 using Base::matrix;
153 using Base::m_error;
154 using Base::m_iterations;
155 using Base::m_info;
156 using Base::m_isInitialized;
157public:
158 typedef _MatrixType MatrixType;
159 typedef typename MatrixType::Scalar Scalar;
160 typedef typename MatrixType::RealScalar RealScalar;
161 typedef _Preconditioner Preconditioner;
162
163public:
164
167
178 template<typename MatrixDerived>
180
182
184 template<typename Rhs,typename Dest>
185 void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
186 {
187 m_iterations = Base::maxIterations();
188 m_error = Base::m_tolerance;
189
190 internal::least_square_conjugate_gradient(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error);
191 m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
192 }
193
194};
195
196} // end namespace Eigen
197
198#endif // EIGEN_LEAST_SQUARE_CONJUGATE_GRADIENT_H
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:144
Index maxIterations() const
Definition: IterativeSolverBase.h:281
A conjugate gradient solver for sparse (or dense) least-square problems.
Definition: LeastSquareConjugateGradient.h:150
LeastSquaresConjugateGradient(const EigenBase< MatrixDerived > &A)
Definition: LeastSquareConjugateGradient.h:179
LeastSquaresConjugateGradient()
Definition: LeastSquareConjugateGradient.h:166
@ Success
Definition: Constants.h:442
@ NoConvergence
Definition: Constants.h:446
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: EigenBase.h:30