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

from __future__ import division
from .minimizer import Minimizer
Martin Reinecke's avatar
Martin Reinecke committed
21
from ..logger import logger
22
from .line_search_strong_wolfe import LineSearchStrongWolfe
23
import numpy as np
24 25


Martin Reinecke's avatar
Martin Reinecke committed
26 27 28
_default_LS = LineSearchStrongWolfe(c2=0.1, preferred_initial_step_size=1.)


29 30
class Yango(Minimizer):
    """ Nonlinear conjugate gradient using curvature
Reimar H Leike's avatar
Reimar H Leike committed
31 32 33 34 35
    The YANGO (Yet Another Nonlinear conjugate Gradient Optimizer)
    uses the curvature to make estimates about suitable descent
    directions. It takes the step that lets it go directly to
    the second order minimum in the subspace spanned by the last
    descent direction and the new gradient.
36 37 38 39 40 41 42 43 44 45 46 47 48 49

    Parameters
    ----------
    controller : IterationController
        Object that decides when to terminate the minimization.

    Notes
    -----
    No restarting procedure has been implemented yet.

    References
    ----------
    """

Martin Reinecke's avatar
Martin Reinecke committed
50
    def __init__(self, controller, line_searcher=_default_LS):
51 52 53 54 55 56 57 58 59 60 61 62
        self._controller = controller
        self._line_searcher = line_searcher

    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
        A_k = energy.curvature
Martin Reinecke's avatar
Martin Reinecke committed
63 64
        energy, success = self._line_searcher.perform_line_search(
            energy, p.vdot(p)/(p.vdot(A_k(p)))*p, f_k_minus_1)
65 66 67
        if not success:
            return energy, controller.ERROR
        A_k = energy.curvature
68
        while True:
Reimar H Leike's avatar
Reimar H Leike committed
69
            r = -energy.gradient
70
            f_k = energy.value
71 72 73 74 75 76
            Ar = A_k(r)
            Ap = A_k(p)
            rAr = r.vdot(Ar)
            pAp = p.vdot(Ap)
            pAr = p.vdot(Ar)
            rAp = r.vdot(Ap)
Reimar H Leike's avatar
Reimar H Leike committed
77 78
            rp = r.vdot(p)
            rr = r.vdot(r)
79
            if rr == 0 or rAr == 0:
Martin Reinecke's avatar
Martin Reinecke committed
80 81
                logger.warning(
                    "Warning: gradient norm 0, assuming convergence!")
82
                return energy, controller.CONVERGED
Martin Reinecke's avatar
Martin Reinecke committed
83
            det = pAp*rAr-np.abs(rAp*pAr)
84
            if det < 0:
Martin Reinecke's avatar
Martin Reinecke committed
85 86 87
                logger.error(
                    "Error: negative determinant ({})".format(det))
                return energy, controller.ERROR
88
            if det == 0:
Martin Reinecke's avatar
Martin Reinecke committed
89 90
                # Try 1D Newton Step
                energy, success = self._line_searcher.perform_line_search(
91 92 93
                    energy, rr/rAr*r, f_k_minus_1)
            else:
                a = (rAr*rp - rAp*rr)/det
94
                b = (pAp*rr - pAr*rp)/det
95
                p = a/b*p+r
Martin Reinecke's avatar
Martin Reinecke committed
96
                energy, success = self._line_searcher.perform_line_search(
97
                    energy, p*b, f_k_minus_1)
98 99
            if not success:
                return energy, controller.ERROR
Reimar H Leike's avatar
Reimar H Leike committed
100
            f_k_minus_1 = f_k
101
            status = self._controller.check(energy)
102 103 104
            if status != controller.CONTINUE:
                return energy, status
            A_k = energy.curvature