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
from __future__ import absolute_import, division, print_function
20

21
from ..compat import *
Martin Reinecke's avatar
Martin Reinecke committed
22
from ..logger import logger
23
from .minimizer import Minimizer
24

Martin Reinecke's avatar
Martin Reinecke committed
25

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

            previous_gamma = gamma