CG.h 4.15 KB
 Tobias Lasser committed Oct 29, 2019 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 #pragma once #include "QuadricProblem.h" #include "Solver.h" #include "LinearOperator.h" namespace elsa { /** * \brief Class implementing the conjugate gradient method * * \author Matthias Wieczorek - initial code * \author David Frank - modularization and modernization * \author Nikola Dinev - rewrite, various enhancements * * CG is an iterative method for minimizing quadric functionals \f$\frac{1}{2} x^tAx - x^tb \f$  Jens Petit committed Oct 31, 2019 17 18 19 20 21 22  * with a symmetric positive-definite operator \f$A \f$. Some common optimization problems, * e.g. weighted least squares or Tikhonov-regularized weighted least squares, can be * reformulated as quadric problems, and solved using CG, even with non-spd operators through * the normal equation. Conversion to the normal equation is performed automatically, where * possible. *  Tobias Lasser committed Oct 29, 2019 23 24 25  * Convergence is considered reached when \f$\| Ax - b \| \leq \epsilon \|Ax_0 - b\| \f$ is * satisfied for some small \f$\epsilon > 0\f$. Here \f$x \f$ denotes the solution * obtained in the last step, and \f$x_0 \f$ denotes the initial guess.  Jens Petit committed Oct 31, 2019 26  *  Tobias Lasser committed Oct 29, 2019 27 28 29 30 31  * References: * https://doi.org/10.6028%2Fjres.049.044 * https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf */ template  Jens Petit committed Oct 31, 2019 32 33  class CG : public Solver {  Tobias Lasser committed Oct 29, 2019 34 35  public: /**  Jens Petit committed Oct 31, 2019 36 37 38  * \brief Constructor for CG, accepting an optimization problem and, optionally, a value for * epsilon *  Tobias Lasser committed Oct 29, 2019 39 40  * \param[in] problem the problem that is supposed to be solved * \param[in] epsilon affects the stopping condition  Jens Petit committed Oct 31, 2019 41 42 43 44 45 46 47  * * If the problem is not a QuadricProblem, a conversion will be attempted. Throws if * conversion fails. See QuadricProblem for details on problems that are convertible to * quadric form. */ explicit CG(const Problem& problem, data_t epsilon = std::numeric_limits::epsilon());  Tobias Lasser committed Oct 29, 2019 48 49  /**  Jens Petit committed Oct 31, 2019 50 51 52  * \brief Constructor for preconditioned CG, accepting an optimization problem, the inverse * of the preconditioner, and, optionally, a value for epsilon *  Tobias Lasser committed Oct 29, 2019 53 54 55  * \param[in] problem the problem that is supposed to be solved * \param[in] preconditionerInverse the inverse of the preconditioner * \param[in] epsilon affects the stopping condition  Jens Petit committed Oct 31, 2019 56 57 58 59 60  * * If the problem is not a QuadricProblem, a conversion will be attempted. Throws if * conversion fails. See QuadricProblem for details on problems that are convertible to * quadric form. */  Tobias Lasser committed Oct 29, 2019 61 62 63 64 65 66  CG(const Problem& problem, const LinearOperator& preconditionerInverse, data_t epsilon = std::numeric_limits::epsilon()); /// default destructor ~CG() override = default;  Nikola Dinev committed Dec 06, 2019 67 68 69 70  protected: /// default copy constructor, hidden from non-derived classes to prevent potential slicing CG(const CG&) = default;  Tobias Lasser committed Oct 29, 2019 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87  private: /// the default number of iterations const index_t _defaultIterations{100}; /// the inverse of the preconditioner (if supplied) std::unique_ptr> _preconditionerInverse{}; /// variable affecting the stopping condition data_t _epsilon; /// lift the base class method getCurrentSolution using Solver::getCurrentSolution; /// lift the base class variable _problem using Solver::_problem; /**  Jens Petit committed Oct 31, 2019 88 89 90 91 92 93 94  * \brief Solve the optimization problem, i.e. apply iterations number of iterations of CG * * \param[in] iterations number of iterations to execute (the default 0 value executes * _defaultIterations of iterations) * * \returns a reference to the current solution */  Jens Petit committed Oct 31, 2019 95  DataContainer& solveImpl(index_t iterations) override;  Tobias Lasser committed Oct 29, 2019 96 97 98 99 100 101 102 103  /// implement the polymorphic clone operation CG* cloneImpl() const override; /// implement the polymorphic comparison operation bool isEqual(const Solver& other) const override; }; } // namespace elsa