Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include <memory>
       4             : 
       5             : #include "DataContainer.h"
       6             : #include "LinearOperator.h"
       7             : #include "Solver.h"
       8             : #include "elsaDefines.h"
       9             : 
      10             : namespace elsa
      11             : {
      12             :     /// @brief Conjugate Gradient for Least Squares Problems
      13             :     ///
      14             :     /// CGLS minimizes:
      15             :     /// \f$ \frac{1}{2} || Ax - b ||_2^2 + \eps^2 || x || \f$
      16             :     /// where \f$A\f$ is an operator, it does not need to be square, symmetric, or positive
      17             :     /// definite. \f$b\f$ is the measured quantity, and \f$\eps\f$ a dampening factors.
      18             :     ///
      19             :     /// If the dampening factor \f$\eps\f$ is non zero, the problem solves a Tikhonov
      20             :     /// problem.
      21             :     ///
      22             :     /// CGLS is equivalent to apply CG to the normal equations \f$A^TAx = A^Tb\f$.
      23             :     /// However, it doesn not need to actually form the normal equation. Primarily,
      24             :     /// this improves the runtime performance, and it further improves the stability
      25             :     /// of the algorithm.
      26             :     ///
      27             :     /// References:
      28             :     /// - https://web.stanford.edu/group/SOL/software/cgls/
      29             :     ///
      30             :     /// @author David Frank
      31             :     template <typename data_t = real_t>
      32             :     class CGLS : public Solver<data_t>
      33             :     {
      34             :     public:
      35             :         /// Scalar alias
      36             :         using Scalar = typename Solver<data_t>::Scalar;
      37             : 
      38             :         /// @brief Construct the necessary form of CGLS using some linear operator and
      39             :         /// the measured data.
      40             :         ///
      41             :         /// @param A linear operator for the problem
      42             :         /// @param b the measured data
      43             :         /// @param eps the dampening factor
      44             :         /// @param tol stopping tolerance
      45             :         CGLS(const LinearOperator<data_t>& A, const DataContainer<data_t>& b,
      46             :              SelfType_t<data_t> eps = 0, SelfType_t<data_t> tol = 1e-4);
      47             : 
      48             :         /// make copy constructor deletion explicit
      49             :         CGLS(const CGLS<data_t>&) = delete;
      50             : 
      51             :         /// default destructor
      52          18 :         ~CGLS() override = default;
      53             : 
      54             :         DataContainer<data_t> setup(std::optional<DataContainer<data_t>> x) override;
      55             : 
      56             :         DataContainer<data_t> step(DataContainer<data_t> x) override;
      57             : 
      58             :         bool shouldStop() const override;
      59             : 
      60             :         std::string formatHeader() const override;
      61             : 
      62             :         std::string formatStep(const DataContainer<data_t>& x) const override;
      63             : 
      64             :     private:
      65             :         std::unique_ptr<LinearOperator<data_t>> A_;
      66             : 
      67             :         DataContainer<data_t> b_;
      68             : 
      69             :         DataContainer<data_t> r_;
      70             : 
      71             :         DataContainer<data_t> s_;
      72             : 
      73             :         DataContainer<data_t> c_;
      74             : 
      75             :         DataContainer<data_t> q_;
      76             : 
      77             :         data_t k_;
      78             : 
      79             :         data_t kold_;
      80             : 
      81             :         data_t damp_{0};
      82             : 
      83             :         data_t tol_{0.0001};
      84             : 
      85             :         /// implement the polymorphic clone operation
      86             :         CGLS<data_t>* cloneImpl() const override;
      87             : 
      88             :         /// implement the polymorphic comparison operation
      89             :         bool isEqual(const Solver<data_t>& other) const override;
      90             :     };
      91             : } // namespace elsa