conjugate_gradient.py 3.68 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
# 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/>.
Theo Steininger's avatar
Theo Steininger committed
13
#
Martin Reinecke's avatar
Martin Reinecke committed
14
# Copyright(C) 2013-2018 Max-Planck-Society
Theo Steininger's avatar
Theo Steininger committed
15
16
17
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
18

19
20
from __future__ import absolute_import, division, print_function
from ..compat import *
Martin Reinecke's avatar
Martin Reinecke committed
21
from .minimizer import Minimizer
Martin Reinecke's avatar
Martin Reinecke committed
22
from ..logger import logger
23

Martin Reinecke's avatar
Martin Reinecke committed
24

Martin Reinecke's avatar
Martin Reinecke committed
25
class ConjugateGradient(Minimizer):
26
27
    """ Implementation of the Conjugate Gradient scheme.

28
29
    It is an iterative method for solving a linear system of equations:
                                    Ax = b
30

31
32
    Parameters
    ----------
Philipp Arras's avatar
Philipp Arras committed
33
    controller : :py:class:`nifty5.IterationController`
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
34
        Object that decides when to terminate the minimization.
35

36
37
    References
    ----------
38
    Jorge Nocedal & Stephen Wright, "Numerical Optimization", Second Edition,
39
40
41
    2006, Springer-Verlag New York
    """

42
    def __init__(self, controller):
Martin Reinecke's avatar
Martin Reinecke committed
43
        self._controller = controller
44

45
    def __call__(self, energy, preconditioner=None):
46
        """ Runs the conjugate gradient minimization.
47
48
49

        Parameters
        ----------
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
50
        energy : Energy object at the starting point of the iteration.
Martin Reinecke's avatar
Martin Reinecke committed
51
            Its metric operator must be independent of position, otherwise
Martin Reinecke's avatar
Martin Reinecke committed
52
            linear conjugate gradient minimization will fail.
53
54
55
        preconditioner : Operator *optional*
            This operator can be provided which transforms the variables of the
            system to improve the conditioning (default: None).
56
57
58

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
59
        QuadraticEnergy
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
60
            state at last point of the iteration
Martin Reinecke's avatar
Martin Reinecke committed
61
        int
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
62
            Can be controller.CONVERGED or controller.ERROR
63
        """
Martin Reinecke's avatar
Martin Reinecke committed
64
        controller = self._controller
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
65
        status = controller.start(energy)
Martin Reinecke's avatar
Martin Reinecke committed
66
        if status != controller.CONTINUE:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
67
            return energy, status
Martin Reinecke's avatar
Martin Reinecke committed
68

69
        r = energy.gradient
70
        d = r if preconditioner is None else preconditioner(r)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
71

Martin Reinecke's avatar
Martin Reinecke committed
72
        previous_gamma = r.vdot(d).real
73
        if previous_gamma == 0:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
74
            return energy, controller.CONVERGED
75

76
        while True:
Martin Reinecke's avatar
Martin Reinecke committed
77
            q = energy.metric(d)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
78
            ddotq = d.vdot(q).real
Martin Reinecke's avatar
Martin Reinecke committed
79
            if ddotq == 0.:
Martin Reinecke's avatar
Martin Reinecke committed
80
                logger.error("Error: ConjugateGradient: ddotq==0.")
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
81
82
                return energy, controller.ERROR
            alpha = previous_gamma/ddotq
83

Martin Reinecke's avatar
Martin Reinecke committed
84
            if alpha < 0:
Martin Reinecke's avatar
Martin Reinecke committed
85
                logger.error("Error: ConjugateGradient: alpha<0.")
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
86
87
                return energy, controller.ERROR

88
            r = r - q*alpha
Martin Reinecke's avatar
Martin Reinecke committed
89

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
90
            energy = energy.at_with_grad(energy.position - alpha*d, r)
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
91

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
92
            s = r if preconditioner is None else preconditioner(r)
93

Martin Reinecke's avatar
Martin Reinecke committed
94
            gamma = r.vdot(s).real
Martin Reinecke's avatar
Martin Reinecke committed
95
            if gamma < 0:
Martin Reinecke's avatar
Martin Reinecke committed
96
                logger.error(
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
97
                    "Positive definiteness of preconditioner violated!")
98
                return energy, controller.ERROR
99
            if gamma == 0:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
100
                return energy, controller.CONVERGED
101

Martin Reinecke's avatar
Martin Reinecke committed
102
            status = controller.check(energy)
103
104
105
            if status != controller.CONTINUE:
                return energy, status

106
            d = d * max(0, gamma/previous_gamma) + s
107
108

            previous_gamma = gamma