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 division
Martin Reinecke's avatar
Martin Reinecke committed
20
from .minimizer import Minimizer
Martin Reinecke's avatar
Martin Reinecke committed
21
22
from ..field import Field
from .. import dobj
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
    ----------
Martin Reinecke's avatar
Martin Reinecke committed
33
    controller : :py:class:`nifty4.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
51
        energy : Energy object at the starting point of the iteration.
            Its curvature 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
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
70
71
        d = r.copy() if preconditioner is None else preconditioner(r)

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
cleanup    
Martin Reinecke committed
77
78
            q = energy.curvature(d)
            ddotq = d.vdot(q).real
Martin Reinecke's avatar
Martin Reinecke committed
79
            if ddotq == 0.:
80
                dobj.mprint("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:
85
                dobj.mprint("Error: ConjugateGradient: alpha<0.")
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
86
87
                return energy, controller.ERROR

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
88
89
            q *= -alpha
            r += q
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:
97
                dobj.mprint(
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

103
            status = self._controller.check(energy)
104
105
106
            if status != controller.CONTINUE:
                return energy, status

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

            previous_gamma = gamma