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 AB_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 AB_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             :      * AB_GMRES solves \f$ min_{y \epsilon \mathbb{R}^m} \left\| b - ABy \right\|\f$, \f$ x = By
      28             :      * \f$.
      29             :      *
      30             :      * \f$ A \f$ can be imagined as a discretization of the forward projector, while \f$ B \f$
      31             :      * represents an matched or unmatched back projector or preconditioner/backprojector, the
      32             :      * algorithm can reach semi-convergence when \f$ A^T \neq B \f$ as long as valid models for
      33             :      * the projector / backprojector pair are chosen
      34             :      *
      35             :      * Convergence is considered reached when \f$ ||r|| <= epsilon \f$
      36             :      *
      37             :      * References:
      38             :      * http://epubs.siam.org/doi/10.1137/0907058
      39             :      * http://epubs.siam.org/doi/10.1137/070696313
      40             :      * http://arxiv.org/abs/2110.01481
      41             :      *
      42             :      * @author Daniel Klitzner - initial code
      43             :      */
      44             :     template <typename data_t = real_t>
      45             :     class AB_GMRES : public Solver<data_t>
      46             :     {
      47             :     public:
      48             :         /// Scalar alias
      49             :         using Scalar = typename Solver<data_t>::Scalar;
      50             : 
      51             :         /**
      52             :          * @brief Constructor for AB_GMRES, accepting a forward projector A, *a sinogram b and an
      53             :          * optional value for epsilon, the backprojector B will be the adjoint *of projector A
      54             :          *
      55             :          * @param[in] projector Discretization of a forward projector
      56             :          * @param[in] sinogram sinogram of measured Data when applying phantom to projector
      57             :          * @param[in] epsilon affects the stopping condition
      58             :          */
      59             :         explicit AB_GMRES(const LinearOperator<data_t>& projector,
      60             :                           const DataContainer<data_t>& sinogram,
      61             :                           data_t epsilon = std::numeric_limits<data_t>::epsilon());
      62             : 
      63             :         /**
      64             :          * @brief Constructor for AB_GMRES, accepting a forward projector A, an *unmatched
      65             :          *backprojector B, a sinogram b and an optional value for *epsilon
      66             :          *
      67             :          * @param[in] projector Discretization of a forward projector A
      68             :          * @param[in] backprojector An optional unmatched/matched backprojector B
      69             :          * @param[in] sinogram sinogram "b" of measured Data when applying phantom to projector
      70             :          * @param[in] epsilon affects the stopping condition
      71             :          */
      72             :         AB_GMRES(const LinearOperator<data_t>& projector,
      73             :                  const LinearOperator<data_t>& backprojector, const DataContainer<data_t>& sinogram,
      74             :                  data_t epsilon = std::numeric_limits<data_t>::epsilon());
      75             : 
      76             :         /// make copy constructor deletion explicit
      77             :         AB_GMRES(const AB_GMRES<data_t>&) = delete;
      78             : 
      79             :         /// default destructor
      80           4 :         ~AB_GMRES() override = default;
      81             : 
      82             :         /**
      83             :          * @brief Solves the given Linear System using AB-GMRES as default, i.e. apply i iterations
      84             :          * with r restarts of AB-GMRES
      85             :          *
      86             :          * @param[in] iterations number of iterations to execute
      87             :          * @param[in] restarts number of restarts to execute after iterations
      88             :          * @param[in] x0 optional approximate starting positions (used for preconditioning)
      89             :          * @return an approximated solution
      90             :          */
      91             :         DataContainer<data_t>
      92             :             solveAndRestart(index_t iterations, index_t restarts,
      93             :                             std::optional<DataContainer<data_t>> x0 = std::nullopt);
      94             : 
      95             :         /**
      96             :          * @brief Solves the given Linear System using AB-GMRES as default, i.e. apply i iterations
      97             :          * of AB-GMRES
      98             :          *
      99             :          * @param[in] iterations number of iterations to execute
     100             :          * @param[in] x0 optional approximate starting positions (used for preconditioning)
     101             :          * @return an approximated solution
     102             :          */
     103             :         DataContainer<data_t>
     104             :             solve(index_t iterations,
     105             :                   std::optional<DataContainer<data_t>> x0 = std::nullopt) override;
     106             : 
     107             :     private:
     108             :         /// Projector A
     109             :         std::unique_ptr<LinearOperator<data_t>> _A;
     110             : 
     111             :         /// Unmatched Backprojector B, used for Preconditioning
     112             :         std::unique_ptr<LinearOperator<data_t>> _B;
     113             : 
     114             :         /// sinogram b
     115             :         DataContainer<data_t> _b;
     116             : 
     117             :         /// variable affecting the stopping condition
     118             :         data_t _epsilon;
     119             : 
     120             :         /// implement the polymorphic clone operation
     121             :         AB_GMRES<data_t>* cloneImpl() const override;
     122             : 
     123             :         /// implement the polymorphic comparison operation
     124             :         bool isEqual(const Solver<data_t>& other) const override;
     125             :     };
     126             : } // namespace elsa