nonlinear_cg.py 3.55 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
14
# Copyright(C) 2013-2019 Max-Planck-Society
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
17

Martin Reinecke's avatar
Martin Reinecke committed
18
from .line_search import LineSearch
19
from .minimizer import Minimizer
20
21
22


class NonlinearCG(Minimizer):
23
    """Nonlinear Conjugate Gradient scheme according to Polak-Ribiere.
Martin Reinecke's avatar
Martin Reinecke committed
24
25

    Algorithm 5.4 from Nocedal & Wright.
Martin Reinecke's avatar
Martin Reinecke committed
26

27
28
29
30
    Parameters
    ----------
    controller : IterationController
        Object that decides when to terminate the minimization.
31
    beta_heuristics : str
Martin Reinecke's avatar
Martin Reinecke committed
32
33
34
35
36
        One of 'Polak-Ribiere', 'Fletcher-Reeves', 'Hestenes-Stiefel' or '5.49'

    Notes
    -----
    No restarting procedure has been implemented yet.
37
38
39
40
41
42
43

    References
    ----------
    Jorge Nocedal & Stephen Wright, "Numerical Optimization", Second Edition,
    2006, Springer-Verlag New York
    """

44
    def __init__(self, controller, beta_heuristics='Polak-Ribiere'):
Martin Reinecke's avatar
Martin Reinecke committed
45
        valid_beta_heuristics = ['Polak-Ribiere', 'Fletcher-Reeves',
46
                                 'Hestenes-Stiefel', "5.49"]
47
        if not (beta_heuristics in valid_beta_heuristics):
Martin Reinecke's avatar
Martin Reinecke committed
48
            raise ValueError("beta heuristics must be either 'Polak-Ribiere', "
Martin Reinecke's avatar
Martin Reinecke committed
49
                             "'Fletcher-Reeves', 'Hestenes-Stiefel, or '5.49'")
50
        self._beta_heuristic = beta_heuristics
51
        self._controller = controller
Martin Reinecke's avatar
Martin Reinecke committed
52
        self._line_searcher = LineSearch(c2=0.1)
53
54
55
56
57
58
59
60
61
62
63
64
65

    def __call__(self, energy):
        controller = self._controller
        status = controller.start(energy)
        if status != controller.CONTINUE:
            return energy, status
        f_k_minus_1 = None

        p = -energy.gradient

        while True:
            grad_old = energy.gradient
            f_k = energy.value
Martin Reinecke's avatar
stage 1    
Martin Reinecke committed
66
67
68
69
            energy, success = self._line_searcher.perform_line_search(
                energy, p, f_k_minus_1)
            if not success:
                return energy, controller.ERROR
70
71
72
73
74
            f_k_minus_1 = f_k
            status = self._controller.check(energy)
            if status != controller.CONTINUE:
                return energy, status
            grad_new = energy.gradient
Martin Reinecke's avatar
Martin Reinecke committed
75

76
            if self._beta_heuristic == 'Hestenes-Stiefel':
77
                # Eq. (5.46) in Nocedal & Wright.
Martin Reinecke's avatar
Martin Reinecke committed
78
79
80
                beta = max(0.0, (grad_new.vdot(grad_new-grad_old) /
                                 (grad_new-grad_old).vdot(p)).real)
            elif self._beta_heuristic == 'Polak-Ribiere':
81
                # Eq. (5.44) in Nocedal & Wright. (with (5.45) additionally)
Martin Reinecke's avatar
Martin Reinecke committed
82
83
                beta = max(0.0, (grad_new.vdot(grad_new-grad_old) /
                                 (grad_old.vdot(grad_old))).real)
84
85
            elif self._beta_heuristic == 'Fletcher-Reeves':
                # Eq. (5.41a) in Nocedal & Wright.
86
                beta = (grad_new.vdot(grad_new)/(grad_old.vdot(grad_old))).real
87
88
89
90
            else:
                # Eq. (5.49) in Nocedal & Wright.
                beta = (grad_new.vdot(grad_new) /
                        ((grad_new-grad_old).vdot(p))).real
91
            p = beta*p - grad_new