Line data    Source code

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