CG.h 4.15 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#pragma once

#include "QuadricProblem.h"
#include "Solver.h"
#include "LinearOperator.h"

namespace elsa
{
    /**
     * \brief Class implementing the conjugate gradient method
     *
     * \author Matthias Wieczorek - initial code
     * \author David Frank - modularization and modernization
     * \author Nikola Dinev - rewrite, various enhancements
     *
     * CG is an iterative method for minimizing quadric functionals \f$ \frac{1}{2} x^tAx - x^tb \f$
Jens Petit's avatar
Jens Petit committed
17
18
19
20
21
22
     * with a symmetric positive-definite operator \f$ A \f$. Some common optimization problems,
     * e.g. weighted least squares or Tikhonov-regularized weighted least squares, can be
     * reformulated as quadric problems, and solved using CG, even with non-spd operators through
     * the normal equation. Conversion to the normal equation is performed automatically, where
     * possible.
     *
23
24
25
     * Convergence is considered reached when \f$ \| Ax - b \| \leq \epsilon \|Ax_0 - b\| \f$ is
     * satisfied for some small \f$ \epsilon > 0\f$. Here \f$ x \f$ denotes the solution
     * obtained in the last step, and \f$ x_0 \f$ denotes the initial guess.
Jens Petit's avatar
Jens Petit committed
26
     *
27
28
29
30
31
     * References:
     * https://doi.org/10.6028%2Fjres.049.044
     * https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf
     */
    template <typename data_t = real_t>
Jens Petit's avatar
Jens Petit committed
32
33
    class CG : public Solver<data_t>
    {
34
35
    public:
        /**
Jens Petit's avatar
Jens Petit committed
36
37
38
         * \brief Constructor for CG, accepting an optimization problem and, optionally, a value for
         * epsilon
         *
39
40
         * \param[in] problem the problem that is supposed to be solved
         * \param[in] epsilon affects the stopping condition
Jens Petit's avatar
Jens Petit committed
41
42
43
44
45
46
47
         *
         * If the problem is not a QuadricProblem, a conversion will be attempted. Throws if
         * conversion fails. See QuadricProblem for details on problems that are convertible to
         * quadric form.
         */
        explicit CG(const Problem<data_t>& problem,
                    data_t epsilon = std::numeric_limits<data_t>::epsilon());
48
49

        /**
Jens Petit's avatar
Jens Petit committed
50
51
52
         * \brief Constructor for preconditioned CG, accepting an optimization problem, the inverse
         * of the preconditioner, and, optionally, a value for epsilon
         *
53
54
55
         * \param[in] problem the problem that is supposed to be solved
         * \param[in] preconditionerInverse the inverse of the preconditioner
         * \param[in] epsilon affects the stopping condition
Jens Petit's avatar
Jens Petit committed
56
57
58
59
60
         *
         * If the problem is not a QuadricProblem, a conversion will be attempted. Throws if
         * conversion fails. See QuadricProblem for details on problems that are convertible to
         * quadric form.
         */
61
62
63
64
65
66
        CG(const Problem<data_t>& problem, const LinearOperator<data_t>& preconditionerInverse,
           data_t epsilon = std::numeric_limits<data_t>::epsilon());

        /// default destructor
        ~CG() override = default;

67
68
69
70
    protected:
        /// default copy constructor, hidden from non-derived classes to prevent potential slicing
        CG(const CG<data_t>&) = default;

71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
    private:
        /// the default number of iterations
        const index_t _defaultIterations{100};

        /// the inverse of the preconditioner (if supplied)
        std::unique_ptr<LinearOperator<data_t>> _preconditionerInverse{};

        /// variable affecting the stopping condition
        data_t _epsilon;

        /// lift the base class method getCurrentSolution
        using Solver<data_t>::getCurrentSolution;

        /// lift the base class variable _problem
        using Solver<data_t>::_problem;

        /**
Jens Petit's avatar
Jens Petit committed
88
89
90
91
92
93
94
         * \brief Solve the optimization problem, i.e. apply iterations number of iterations of CG
         *
         * \param[in] iterations number of iterations to execute (the default 0 value executes
         * _defaultIterations of iterations)
         *
         * \returns a reference to the current solution
         */
95
        DataContainer<data_t>& solveImpl(index_t iterations) override;
96
97
98
99
100
101
102
103

        /// implement the polymorphic clone operation
        CG<data_t>* cloneImpl() const override;

        /// implement the polymorphic comparison operation
        bool isEqual(const Solver<data_t>& other) const override;
    };
} // namespace elsa