conjugate_gradient.py 3.66 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 division
Martin Reinecke's avatar
Martin Reinecke committed
20
from .minimizer import Minimizer
Martin Reinecke's avatar
Martin Reinecke committed
21
from ..logger import logger
22

Martin Reinecke's avatar
Martin Reinecke committed
23

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

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

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

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

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

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

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

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

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

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

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

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

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
87
            q *= -alpha
Martin Reinecke's avatar
Martin Reinecke committed
88
            r = r + q
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

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
106
107
            d *= max(0, gamma/previous_gamma)
            d += s
108
109

            previous_gamma = gamma