Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include <memory>
       4             : #include <optional>
       5             : 
       6             : #include "Solver.h"
       7             : #include "LinearOperator.h"
       8             : #include <Eigen/Core>
       9             : 
      10             : namespace elsa
      11             : {
      12             :     /**
      13             :      * @brief Class implementing the BA_GMRES method for solving with an unmatched Back Projector
      14             :      *
      15             :      * GMRES is an iterative method for solving an indefinite nonsymmetric linear system \f$ Ax = b
      16             :      * \f$ It solves the system by minimizing a least squares Problem for \f$ J(y) = \left\| \beta
      17             :      * {e_{1}} - \bar{H}y \right\| \f$ where \f$ \beta \f$ is \f$ \left\| f - Ax \right\| \f$, \f$
      18             :      * {e_{1}} \f$ is the first column of the \f$ (k + 1) \times (k + 1) \f$ identity matrix and \f$
      19             :      * \bar{H} \f$ is \f$ (k + 1) \times (k) \f$ matrix whose only nonzero entries are generated by
      20             :      * the algorithm.
      21             :      *
      22             :      * We implement an expanded BA_GMRES method for solving large sparse least square problems. Here
      23             :      * we define \f$ A^T = B \f$ for our matrix \f$ A \epsilon \mathbb{R}^{m \times n} \f$ with
      24             :      * \f$ B \epsilon \mathbb{R}^{n \times m} \f$. We solve the so-called unmatched normal equations
      25             :      * \f$ BAx = Bb \f$ and \f$ ABy = b \f$  where \f$ x = By \f$:
      26             :      *
      27             :      * BA_GMRES solves \f$ min_{x \epsilon \mathbb{R}^m} \left\| Bb - BAx \right\|\f$.
      28             :      *
      29             :      * \f$ A \f$ can be imagined as a discretization of the forward projector, while \f$ B \f$
      30             :      * represents an matched or unmatched back projector or preconditioner/backprojector, the
      31             :      * algorithm can reach semi-convergence when \f$ A^T \neq B \f$ as long as valid models for
      32             :      * the projector / backprojector pair are chosen
      33             :      *
      34             :      * Convergence is considered reached when \f$ ||r|| <= epsilon \f$
      35             :      *
      36             :      * References:
      37             :      * http://epubs.siam.org/doi/10.1137/0907058
      38             :      * http://epubs.siam.org/doi/10.1137/070696313
      39             :      * http://arxiv.org/abs/2110.01481
      40             :      *
      41             :      * @author Daniel Klitzner - initial code
      42             :      */
      43             :     template <typename data_t = real_t>
      44             :     class BA_GMRES : public Solver<data_t>
      45             :     {
      46             :     public:
      47             :         /// Scalar alias
      48             :         using Scalar = typename Solver<data_t>::Scalar;
      49             : 
      50             :         /**
      51             :          * @brief Constructor for BA_GMRES, accepting a forward projector A, *a sinogram b and an
      52             :          *optional value for epsilon, the backprojector B will be the adjoint *of projector A
      53             :          *
      54             :          * @param[in] projector Discretization of a forward projector
      55             :          * @param[in] sinogram sinogram of measured Data when applying phantom to projector
      56             :          * @param[in] epsilon affects the stopping condition
      57             :          **/
      58             :         explicit BA_GMRES(const LinearOperator<data_t>& projector,
      59             :                           const DataContainer<data_t>& sinogram,
      60             :                           data_t epsilon = std::numeric_limits<data_t>::epsilon());
      61             : 
      62             :         /**
      63             :          * @brief Constructor for BA_GMRES, accepting a forward projector A, an *unmatched
      64             :          *backprojector B, a sinogram b and an optional value for *epsilon
      65             :          *
      66             :          * @param[in] projector Discretization of a forward projector A
      67             :          * @param[in] backprojector An optional unmatched/matched backprojector B
      68             :          * @param[in] sinogram sinogram "b" of measured Data when applying phantom to projector
      69             :          * @param[in] epsilon affects the stopping condition
      70             :          **/
      71             :         BA_GMRES(const LinearOperator<data_t>& projector,
      72             :                  const LinearOperator<data_t>& backprojector, const DataContainer<data_t>& sinogram,
      73             :                  data_t epsilon = std::numeric_limits<data_t>::epsilon());
      74             : 
      75             :         /// make copy constructor deletion explicit
      76             :         BA_GMRES(const BA_GMRES<data_t>&) = delete;
      77             : 
      78             :         /// default destructor
      79           4 :         ~BA_GMRES() override = default;
      80             : 
      81             :         /**
      82             :          * @brief Solves the given Linear System using BA-GMRES as default, i.e. apply i iterations
      83             :          *with r restarts of BA-GMRES
      84             :          *
      85             :          * @param[in] iterations number of iterations to execute
      86             :          * @param[in] restarts number of restarts to execute after iterations
      87             :          * @param[in] x0 optional approximate starting positions (used for preconditioning)
      88             :          * @return an approximated solution
      89             :          **/
      90             :         DataContainer<data_t>
      91             :             solveAndRestart(index_t iterations, index_t restarts,
      92             :                             std::optional<DataContainer<data_t>> x0 = std::nullopt);
      93             : 
      94             :         /**
      95             :          * @brief Solves the given Linear System using BA-GMRES as default, i.e. apply i iterations
      96             :          *of BA-GMRES
      97             :          *
      98             :          * @param[in] iterations number of iterations to execute
      99             :          * @param[in] x0 optional approximate starting positions (used for preconditioning)
     100             :          * @return an approximated solution
     101             :          **/
     102             :         DataContainer<data_t>
     103             :             solve(index_t iterations,
     104             :                   std::optional<DataContainer<data_t>> x0 = std::nullopt) override;
     105             : 
     106             :     private:
     107             :         /// Projector A
     108             :         std::unique_ptr<LinearOperator<data_t>> _A;
     109             : 
     110             :         /// Unmatched Backprojector B, used for Preconditioning
     111             :         std::unique_ptr<LinearOperator<data_t>> _B;
     112             : 
     113             :         /// sinogram b
     114             :         DataContainer<data_t> _b;
     115             : 
     116             :         /// variable affecting the stopping condition
     117             :         data_t _epsilon;
     118             : 
     119             :         /// implement the polymorphic clone operation
     120             :         BA_GMRES<data_t>* cloneImpl() const override;
     121             : 
     122             :         /// implement the polymorphic comparison operation
     123             :         bool isEqual(const Solver<data_t>& other) const override;
     124             :     };
     125             : } // namespace elsa