Line data    Source code

       1             : 
       2             : #include <optional>
       3             : 
       4             : #include "Functional.h"
       5             : #include "Solver.h"
       6             : #include "LinearOperator.h"
       7             : #include "StrongTypes.h"
       8             : #include "MaybeUninitialized.hpp"
       9             : #include "ProximalOperator.h"
      10             : #include "LineSearchMethod.h"
      11             : #include "FixedStepSize.h"
      12             : 
      13             : namespace elsa
      14             : {
      15             :     /**
      16             :      * @brief Proximal Optimized Gradient Method (POGM)
      17             :      *
      18             :      * POGM minimizes function of the the same for as PGD. See the documentation there.
      19             :      *
      20             :      * The update steps for POGM are:
      21             :      * @f[
      22             :      *    \begin{split}
      23             :      *        \theta_k &= \begin{cases}%
      24             :      *                \frac{1}{2}(1 + \sqrt{4 \theta_{k-1}^2 + 1}) & \text{ if } 2 \leq k < N \\%
      25             :      *                \frac{1}{2}(1 + \sqrt{8 \theta_{k-1}^2 + 1}) & \text{ if } k = N        %
      26             :      *        \end{cases} \\%
      27             :      *        \gamma_k &= \frac{1}{L} \frac{2\theta_{k-1} + \theta_k - 1}{\theta_k} \\%
      28             :      *        \omega_k &= x_{k+1} - \frac{1}{L} \nabla f(x_{k-1}) \\%
      29             :      *        z_k &= \omega_k +%
      30             :      *            \underset{Nesterov}{\underbrace{\frac{\theta_{k-1} - 1}{\theta_k}(\omega_k +
      31             :      *\omega_{k-1})}}+% \underset{OGM}{\underbrace{\frac{\theta_{k-1}}{\theta_k}(\omega_k -
      32             :      *x_{k-1})}}+% \underset{POGM}{\underbrace{\frac{\theta_{k-1} - 1}{L
      33             :      *\gamma_{k-1}\theta_k}(z_{k-1} - x_{k-1})}}\\% x_k &= \operatorname{prox}_{\gamma_k g}(z_k)
      34             :      * \end{split}
      35             :      * @f]
      36             :      *
      37             :      * References:
      38             :      * - https://arxiv.org/abs/1512.07516
      39             :      *
      40             :      * @see For a more detailed discussion of the type of problem for this solver,
      41             :      * see PGD.
      42             :      *
      43             :      * @author
      44             :      * - David Frank - Initial implementation
      45             :      *
      46             :      * @tparam data_t data type for the domain and range of the problem, defaulting to real_t
      47             :      */
      48             :     template <typename data_t = real_t>
      49             :     class POGM : public Solver<data_t>
      50             :     {
      51             :     public:
      52             :         /// Scalar alias
      53             :         using Scalar = typename Solver<data_t>::Scalar;
      54             : 
      55             :         /// Construct POGM with a least squares data fidelity term
      56             :         ///
      57             :         /// @note The step length for least squares can be chosen to be dependend on the Lipschitz
      58             :         /// constant. Compute it using `powerIterations(adjoint(A) * A)`. Depending on which
      59             :         /// literature, both \f$ \frac{2}{L} \f$ and \f$ \frac{1}{L}\f$. If mu is not given, the
      60             :         /// step length is chosen by default, the computation of the power method might be
      61             :         /// expensive.
      62             :         ///
      63             :         /// @param A the operator for the least squares data term
      64             :         /// @param b the measured data for the least squares data term
      65             :         /// @param g prox-friendly function
      66             :         /// @param mu the step length
      67             :         POGM(const LinearOperator<data_t>& A, const DataContainer<data_t>& b,
      68             :              const Functional<data_t>& h, std::optional<data_t> mu = std::nullopt,
      69             :              data_t epsilon = std::numeric_limits<data_t>::epsilon());
      70             : 
      71             :         /// Construct POGM with a weighted least squares data fidelity term
      72             :         ///
      73             :         /// @note The step length for least squares can be chosen to be dependend on the Lipschitz
      74             :         /// constant. Compute it using `powerIterations(adjoint(A) * A)`. Depending on which
      75             :         /// literature, both \f$ \frac{2}{L} \f$ and \f$ \frac{1}{L}\f$. If mu is not given, the
      76             :         /// step length is chosen by default, the computation of the power method might be
      77             :         /// expensive.
      78             :         ///
      79             :         /// @param A the operator for the least squares data term
      80             :         /// @param b the measured data for the least squares data term
      81             :         /// @param W the weights (usually `counts / I0`)
      82             :         /// @param prox the proximal operator for g
      83             :         /// @param mu the step length
      84             :         POGM(const LinearOperator<data_t>& A, const DataContainer<data_t>& b,
      85             :              const DataContainer<data_t>& W, const Functional<data_t>& h,
      86             :              std::optional<data_t> mu = std::nullopt,
      87             :              data_t epsilon = std::numeric_limits<data_t>::epsilon());
      88             : 
      89             :         /// Construct POGM with a given data fidelity term
      90             :         POGM(const Functional<data_t>& g, const Functional<data_t>& h, data_t mu,
      91             :              data_t epsilon = std::numeric_limits<data_t>::epsilon());
      92             : 
      93             :         /// Construct POGM with a least squares data fidelity term
      94             :         ///
      95             :         /// @note The step length for least squares can be chosen to be dependend on the Lipschitz
      96             :         /// constant. Compute it using `powerIterations(adjoint(A) * A)`. Depending on which
      97             :         /// literature, both \f$ \frac{2}{L} \f$ and \f$ \frac{1}{L}\f$.
      98             :         ///
      99             :         /// @param A the operator for the least squares data term
     100             :         /// @param b the measured data for the least squares data term
     101             :         /// @param g prox-friendly function
     102             :         /// @param lineSearchMethod the line search method to find the step size at
     103             :         /// each iteration
     104             :         POGM(const LinearOperator<data_t>& A, const DataContainer<data_t>& b,
     105             :              const Functional<data_t>& h, const LineSearchMethod<data_t>& lineSearchMethod,
     106             :              data_t epsilon = std::numeric_limits<data_t>::epsilon());
     107             : 
     108             :         /// Construct POGM with a weighted least squares data fidelity term
     109             :         ///
     110             :         /// @note The step length for least squares can be chosen to be dependend on the Lipschitz
     111             :         /// constant. Compute it using `powerIterations(adjoint(A) * A)`. Depending on which
     112             :         /// literature, both \f$ \frac{2}{L} \f$ and \f$ \frac{1}{L}\f$.
     113             :         ///
     114             :         /// @param A the operator for the least squares data term
     115             :         /// @param b the measured data for the least squares data term
     116             :         /// @param W the weights (usually `counts / I0`)
     117             :         /// @param prox the proximal operator for g
     118             :         /// @param lineSearchMethod the line search method to find the step size at
     119             :         /// each iteration
     120             :         POGM(const LinearOperator<data_t>& A, const DataContainer<data_t>& b,
     121             :              const DataContainer<data_t>& W, const Functional<data_t>& h,
     122             :              const LineSearchMethod<data_t>& lineSearchMethod,
     123             :              data_t epsilon = std::numeric_limits<data_t>::epsilon());
     124             : 
     125             :         /// Construct POGM with a given data fidelity term
     126             :         POGM(const Functional<data_t>& g, const Functional<data_t>& h,
     127             :              const LineSearchMethod<data_t>& lineSearchMethod,
     128             :              data_t epsilon = std::numeric_limits<data_t>::epsilon());
     129             : 
     130             :         /// default destructor
     131           8 :         ~POGM() override = default;
     132             : 
     133             :         /**
     134             :          * @brief Solve the optimization problem, i.e. apply iterations number of iterations of
     135             :          * POGM
     136             :          *
     137             :          * @param[in] iterations number of iterations to execute
     138             :          *
     139             :          * @returns the approximated solution
     140             :          */
     141             :         DataContainer<data_t>
     142             :             solve(index_t iterations,
     143             :                   std::optional<DataContainer<data_t>> x0 = std::nullopt) override;
     144             : 
     145             :     protected:
     146             :         /// implement the polymorphic clone operation
     147             :         auto cloneImpl() const -> POGM<data_t>* override;
     148             : 
     149             :         /// implement the polymorphic comparison operation
     150             :         auto isEqual(const Solver<data_t>& other) const -> bool override;
     151             : 
     152             :     private:
     153             :         /// The LASSO optimization problem
     154             :         std::unique_ptr<Functional<data_t>> g_;
     155             : 
     156             :         std::unique_ptr<Functional<data_t>> h_;
     157             : 
     158             :         // ProximalOperator<data_t> prox_;
     159             : 
     160             :         /// the line search method
     161             :         std::unique_ptr<LineSearchMethod<data_t>> lineSearchMethod_;
     162             : 
     163             :         /// variable affecting the stopping condition
     164             :         data_t epsilon_;
     165             :     };
     166             : } // namespace elsa