merge better_minimizers

nifty/energies/__init__.py

@@ -17,5 +17,6 @@
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
 from .energy import Energy
+from .quadratic_energy import QuadraticEnergy
 from .line_energy import LineEnergy
 from .memoization import memo

nifty/energies/energy.py

@@ -17,6 +17,7 @@
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
 from ..nifty_meta import NiftyMeta
+from .memoization import memo
 
 from keepers import Loggable
 from future.utils import with_metaclass
@@ -41,7 +42,7 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
     value : np.float
         The value of the energy functional at given `position`.
     gradient : Field
-        The gradient at given `position` in parameter direction.
+        The gradient at given `position`.
     curvature : LinearOperator, callable
         A positive semi-definite operator or function describing the curvature
         of the potential at the given `position`.
@@ -108,12 +109,32 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
     @property
     def gradient(self):
         """
-        The gradient at given `position` in parameter direction.
+        The gradient at given `position`.
 
         """
 
         raise NotImplementedError
 
+    @property
+    @memo
+    def gradient_norm(self):
+        """
+        The length of the gradient at given `position`.
+
+        """
+
+        return self.gradient.norm()
+
+    @property
+    @memo
+    def gradient_infnorm(self):
+        """
+        The infinity norm of the gradient at given `position`.
+
+        """
+
+        return abs(self.gradient).max()
+
     @property
     def curvature(self):
         """

nifty/energies/quadratic_energy.py

+from .energy import Energy
+from .memoization import memo
+
+
+class QuadraticEnergy(Energy):
+    """The Energy for a quadratic form.
+    The most important aspect of this energy is that its curvature must be
+    position-independent.
+    """
+
+    def __init__(self, position, A, b, grad=None):
+        super(QuadraticEnergy, self).__init__(position=position)
+        self._A = A
+        self._b = b
+        if grad is not None:
+            self._Ax = grad + self._b
+        else:
+            self._Ax = self._A(self.position)
+
+    def at(self, position):
+        return self.__class__(position=position, A=self._A, b=self._b)
+
+    def at_with_grad(self, position, grad):
+        return self.__class__(position=position, A=self._A, b=self._b,
+                              grad=grad)
+
+    @property
+    @memo
+    def value(self):
+        return 0.5*self.position.vdot(self._Ax) - self._b.vdot(self.position)
+
+    @property
+    @memo
+    def gradient(self):
+        return self._Ax - self._b
+
+    @property
+    def curvature(self):
+        return self._A

nifty/minimization/__init__.py

@@ -17,7 +17,11 @@
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
 from .line_searching import *
+from .iteration_controller import IterationController
+from .default_iteration_controller import DefaultIterationController
+from .minimizer import Minimizer
 from .conjugate_gradient import ConjugateGradient
+from .nonlinear_cg import NonlinearCG
 from .descent_minimizer import DescentMinimizer
 from .steepest_descent import SteepestDescent
 from .vl_bfgs import VL_BFGS

nifty/minimization/conjugate_gradient.py

@@ -19,10 +19,10 @@
 from __future__ import division
 import numpy as np
 
-from keepers import Loggable
+from .minimizer import Minimizer
 
 
-class ConjugateGradient(Loggable, object):
+class ConjugateGradient(Minimizer):
     """ Implementation of the Conjugate Gradient scheme.
 
     It is an iterative method for solving a linear system of equations:
@@ -30,43 +30,11 @@ class ConjugateGradient(Loggable, object):
 
     Parameters
     ----------
-    convergence_tolerance : float *optional*
-        Tolerance specifying the case of convergence. (default: 1E-4)
-    convergence_level : integer *optional*
-        Number of times the tolerance must be undershot before convergence
-        is reached. (default: 3)
-    iteration_limit : integer *optional*
-        Maximum number of iterations performed (default: None).
-    reset_count : integer *optional*
-        Number of iterations after which to restart; i.e., forget previous
-        conjugated directions (default: None).
+    controller : IterationController
+        Object that decides when to terminate the minimization.
     preconditioner : Operator *optional*
         This operator can be provided which transforms the variables of the
         system to improve the conditioning (default: None).
-    callback : callable *optional*
-        Function f(energy, iteration_number) supplied by the user to perform
-        in-situ analysis at every iteration step. When being called the
-        current energy and iteration_number are passed. (default: None)
-
-    Attributes
-    ----------
-    convergence_tolerance : float
-        Tolerance specifying the case of convergence.
-    convergence_level : integer
-        Number of times the tolerance must be undershot before convergence
-        is reached. (default: 3)
-    iteration_limit : integer
-        Maximum number of iterations performed.
-    reset_count : integer
-        Number of iterations after which to restart; i.e., forget previous
-        conjugated directions.
-    preconditioner : function
-        This operator can be provided which transforms the variables of the
-        system to improve the conditioning (default: None).
-    callback : callable
-        Function f(energy, iteration_number) supplied by the user to perform
-        in-situ analysis at every iteration step. When being called the
-        current energy and iteration_number are passed. (default: None)
 
     References
     ----------
@@ -75,140 +43,73 @@ class ConjugateGradient(Loggable, object):
 
     """
 
-    def __init__(self, convergence_tolerance=1E-4, convergence_level=3,
-                 iteration_limit=None, reset_count=None,
-                 preconditioner=None, callback=None):
-
-        self.convergence_tolerance = np.float(convergence_tolerance)
-        self.convergence_level = np.float(convergence_level)
-
-        if iteration_limit is not None:
-            iteration_limit = int(iteration_limit)
-        self.iteration_limit = iteration_limit
+    def __init__(self, controller, preconditioner=None):
+        self._preconditioner = preconditioner
+        self._controller = controller
 
-        if reset_count is not None:
-            reset_count = int(reset_count)
-        self.reset_count = reset_count
-
-        if preconditioner is None:
-            preconditioner = lambda z: z
-
-        self.preconditioner = preconditioner
-        self.callback = callback
-
-    def __call__(self, A, b, x0):
+    def __call__(self, energy):
         """ Runs the conjugate gradient minimization.
-        For `Ax = b` the variable `x` is infered.
 
         Parameters
         ----------
-        A : Operator
-            Operator `A` applicable to a Field.
-        b : Field
-            Result of the operation `A(x)`.
-        x0 : Field
-            Starting guess for the minimization.
+        energy : Energy object at the starting point of the iteration.
+            Its curvature operator must be independent of position, otherwise
+            linear conjugate gradient minimization will fail.
 
         Returns
         -------
-        x : Field
-            Latest `x` of the minimization.
-        convergence : integer
-            Latest convergence level indicating whether the minimization
-            has converged or not.
+        energy : QuadraticEnergy
+            state at last point of the iteration
+        status : integer
+            Can be controller.CONVERGED or controller.ERROR
 
         """
 
-        r = b - A(x0)
-        d = self.preconditioner(r)
+        controller = self._controller
+        status = controller.start(energy)
+        if status != controller.CONTINUE:
+            return energy, status
+
+        r = -energy.gradient
+        if self._preconditioner is not None:
+            d = self._preconditioner(r)
+        else:
+            d = r.copy()
         previous_gamma = (r.vdot(d)).real
         if previous_gamma == 0:
-            self.logger.info("The starting guess is already perfect solution "
-                             "for the inverse problem.")
-            return x0, self.convergence_level+1
-        norm_b = np.sqrt((b.vdot(b)).real)
-        x = x0.copy()
-        convergence = 0
-        iteration_number = 1
-        self.logger.info("Starting conjugate gradient.")
-
-        beta = np.inf
-        delta = np.inf
+            return energy, controller.CONVERGED
 
         while True:
-            if self.callback is not None:
-                self.callback(x, iteration_number)
+            q = energy.curvature(d)
+            ddotq = d.vdot(q).real
+            if ddotq==0.:
+                self.logger.error("Alpha became infinite! Stopping.")
+                return energy, controller.ERROR
+            alpha = previous_gamma/ddotq
 
-            q = A(d)
-            alpha = previous_gamma/d.vdot(q).real
+            if alpha < 0:
+                self.logger.warn("Positive definiteness of A violated!")
+                return energy, controller.ERROR
 
-            if not np.isfinite(alpha):
-                self.logger.error(
-                        "Alpha became infinite! Stopping. Iteration : %08u   "
-                        "alpha = %3.1E   beta = %3.1E   delta = %3.1E" %
-                        (iteration_number, alpha, beta, delta))
-                return x0, 0
+            r -= q * alpha
+            energy = energy.at_with_grad(energy.position+d*alpha,-r)
 
-            x += d * alpha
+            status = self._controller.check(energy)
+            if status != controller.CONTINUE:
+                return energy, status
 
-            reset = False
-            if alpha < 0:
-                self.logger.warn("Positive definiteness of A violated!")
-                reset = True
-            if self.reset_count is not None:
-                reset += (iteration_number % self.reset_count == 0)
-            if reset:
-                self.logger.info("Computing accurate residuum.")
-                r = b - A(x)
+            if self._preconditioner is not None:
+                s = self._preconditioner(r)
             else:
-                r -= q * alpha
+                s = r
 
-            s = self.preconditioner(r)
             gamma = r.vdot(s).real
-
             if gamma < 0:
-                self.logger.warn("Positive definitness of preconditioner "
-                                 "violated!")
-
-            beta = max(0, gamma/previous_gamma)
-
-            delta = np.sqrt(gamma)/norm_b
-
-            self.logger.debug("Iteration : %08u   alpha = %3.1E   "
-                              "beta = %3.1E   delta = %3.1E" %
-                              (iteration_number, alpha, beta, delta))
-
+                self.logger.warn(
+                    "Positive definiteness of preconditioner violated!")
             if gamma == 0:
-                convergence = self.convergence_level+1
-                self.logger.info(
-                        "Reached infinite convergence. Iteration : %08u   "
-                        "alpha = %3.1E   beta = %3.1E   delta = %3.1E" %
-                        (iteration_number, alpha, beta, delta))
-                break
-            elif abs(delta) < self.convergence_tolerance:
-                convergence += 1
-                self.logger.info("Updated convergence level to: %u" %
-                                 convergence)
-                if convergence == self.convergence_level:
-                    self.logger.info(
-                        "Reached target convergence level. Iteration : %08u   "
-                        "alpha = %3.1E   beta = %3.1E   delta = %3.1E" %
-                        (iteration_number, alpha, beta, delta))
-                    break
-            else:
-                convergence = max(0, convergence-1)
+                return energy, controller.CONVERGED
 
-            if self.iteration_limit is not None:
-                if iteration_number == self.iteration_limit:
-                    self.logger.info(
-                        "Reached iteration limit. Iteration : %08u   "
-                        "alpha = %3.1E   beta = %3.1E   delta = %3.1E" %
-                        (iteration_number, alpha, beta, delta))
-                    break
+            d = s + d * max(0, gamma/previous_gamma)
 
-            d = s + d * beta
-
-            iteration_number += 1
             previous_gamma = gamma
-
-        return x, convergence

nifty/minimization/default_iteration_controller.py

+# 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-2017 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 print_function
+from .iteration_controller import IterationController
+
+class DefaultIterationController(IterationController):
+    def __init__ (self, tol_abs_gradnorm=None, tol_rel_gradnorm=None,
+                  convergence_level=1, iteration_limit=None, name=None,
+                  verbose=None):
+        super(DefaultIterationController, self).__init__()
+        self._tol_abs_gradnorm = tol_abs_gradnorm
+        self._tol_rel_gradnorm = tol_rel_gradnorm
+        self._convergence_level = convergence_level
+        self._iteration_limit = iteration_limit
+        self._name = name
+        self._verbose = verbose
+
+    def start(self, energy):
+        self._itcount = -1
+        self._ccount = 0
+        if self._tol_rel_gradnorm is not None:
+            self._tol_rel_gradnorm_now = self._tol_rel_gradnorm \
+                                       * energy.gradient_norm
+        return self.check(energy)
+
+    def check(self, energy):
+        self._itcount += 1
+
+        if self._tol_abs_gradnorm is not None:
+            if energy.gradient_norm <= self._tol_abs_gradnorm:
+                self._ccount += 1
+        if self._tol_rel_gradnorm is not None:
+            if energy.gradient_norm <= self._tol_rel_gradnorm_now:
+                self._ccount += 1
+
+        # report
+        if self._verbose:
+            msg = ""
+            if self._name is not None:
+                msg += self._name+":"
+            msg += " Iteration #" + str(self._itcount)
+            msg += " gradnorm=" + str(energy.gradient_norm)
+            msg += " convergence level=" + str(self._ccount)
+            print (msg)
+            self.logger.info(msg)
+
+        # Are we done?
+        if self._iteration_limit is not None:
+            if self._itcount >= self._iteration_limit:
+                return self.CONVERGED
+        if self._ccount >= self._convergence_level:
+            return self.CONVERGED
+
+        return self.CONTINUE

nifty/minimization/descent_minimizer.py

@@ -18,17 +18,14 @@
 
 from __future__ import division
 import abc
-from ..nifty_meta import NiftyMeta
-
 import numpy as np
 
-from keepers import Loggable
-
+from .minimizer import Minimizer
 from .line_searching import LineSearchStrongWolfe
 from future.utils import with_metaclass
 
 
-class DescentMinimizer(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {}))):
+class DescentMinimizer(Minimizer):
     """ A base class used by gradient methods to find a local minimum.
 
     Descent minimization methods are used to find a local minimum of a scalar
@@ -38,60 +35,18 @@ class DescentMinimizer(with_metaclass(NiftyMeta, type('NewBase', (Loggable, obje
 
     Parameters
     ----------
+    controller : IterationController
+        Object that decides when to terminate the minimization.
     line_searcher : callable *optional*
         Function which infers the step size in the descent direction
         (default : LineSearchStrongWolfe()).
-    callback : callable *optional*
-        Function f(energy, iteration_number) supplied by the user to perform
-        in-situ analysis at every iteration step. When being called the
-        current energy and iteration_number are passed. (default: None)
-    convergence_tolerance : float *optional*
-        Tolerance specifying the case of convergence. (default: 1E-4)
-    convergence_level : integer *optional*
-        Number of times the tolerance must be undershot before convergence
-        is reached. (default: 3)
-    iteration_limit : integer *optional*
-        Maximum number of iterations performed (default: None).
-
-    Attributes
-    ----------
-    convergence_tolerance : float
-        Tolerance specifying the case of convergence.
-    convergence_level : integer
-        Number of times the tolerance must be undershot before convergence
-        is reached. (default: 3)
-    iteration_limit : integer
-        Maximum number of iterations performed.
-    line_searcher : LineSearch
-        Function which infers the optimal step size for functional minization
-        given a descent direction.
-    callback : function
-        Function f(energy, iteration_number) supplied by the user to perform
-        in-situ analysis at every iteration step. When being called the
-        current energy and iteration_number are passed.
-
-    Notes
-    ------
-    The callback function can be used to externally stop the minimization by
-    raising a `StopIteration` exception.
-    Check `get_descent_direction` of a derived class for information on the
-    concrete minization scheme.
 
     """
 
-    def __init__(self, line_searcher=LineSearchStrongWolfe(), callback=None,
-                 convergence_tolerance=1E-4, convergence_level=3,
-                 iteration_limit=None):
-
-        self.convergence_tolerance = np.float(convergence_tolerance)
-        self.convergence_level = np.int(convergence_level)
-
-        if iteration_limit is not None:
-            iteration_limit = int(iteration_limit)
-        self.iteration_limit = iteration_limit
-
+    def __init__(self, controller, line_searcher=LineSearchStrongWolfe()):
+        super(DescentMinimizer, self).__init__()
+        self._controller = controller
         self.line_searcher = line_searcher
-        self.callback = callback
 
     def __call__(self, energy):
         """ Performs the minimization of the provided Energy functional.
@@ -106,43 +61,29 @@ class DescentMinimizer(with_metaclass(NiftyMeta, type('NewBase', (Loggable, obje
         -------
         energy : Energy object
             Latest `energy` of the minimization.
-        convergence : integer
-            Latest convergence level indicating whether the minimization
-            has converged or not.
+        status : integer
+            Can be controller.CONVERGED or controller.ERROR
 
         Note
         ----
         The minimization is stopped if
-            * the callback function raises a `StopIteration` exception,
+            * the controller returns controller.CONVERGED or controller.ERROR,
             * a perfectly flat point is reached,
             * according to the line-search the minimum is found,
-            * the target convergence level is reached,
-            * the iteration limit is reached.
 
         """
 
-        convergence = 0
         f_k_minus_1 = None
-        iteration_number = 1
+        controller = self._controller
+        status = controller.start(energy)
+        if status != controller.CONTINUE:
+            return energy, status
 
         while True:
-            if self.callback is not None:
-                try:
-                    self.callback(energy, iteration_number)
-                except StopIteration:
-                    self.logger.info("Minimization was stopped by callback "
-                                     "function.")
-                    break
-
-            # compute the the gradient for the current location
-            gradient = energy.gradient
-            gradient_norm = gradient.norm()
-
             # check if position is at a flat point
-            if gradient_norm == 0:
+            if energy.gradient_norm == 0:
                 self.logger.info("Reached perfectly flat point. Stopping.")
-                convergence = self.convergence_level+2
-                break
+                return energy, controller.CONVERGED
 
             # current position is encoded in energy object
             descent_direction = self.get_descent_direction(energy)
@@ -157,47 +98,20 @@ class DescentMinimizer(with_metaclass(NiftyMeta, type('NewBase', (Loggable, obje
             except RuntimeError:
                 self.logger.warn(
                         "Stopping because of RuntimeError in line-search")
-                break
+                return energy, controller.ERROR
 
             f_k_minus_1 = energy.value
-            f_k = new_energy.value
-            delta = (abs(f_k-f_k_minus_1) /
-                     max(abs(f_k), abs(f_k_minus_1), 1.))
             # check if new energy value is bigger than old energy value
             if (new_energy.value - energy.value) > 0:
                 self.logger.info("Line search algorithm returned a new energy "
                                  "that was larger than the old one. Stopping.")
-                break
+                return energy, controller.ERROR
 
             energy = new_energy
-            # check convergence
-            self.logger.debug("Iteration:%08u "
-                              "delta=%3.1E energy=%3.1E" %
-                              (iteration_number, delta,
-                               energy.value))
-            if delta == 0:
-                convergence = self.convergence_level + 2
-                self.logger.info("Found minimum according to line-search. "
-                                 "Stopping.")
-                break
-            elif delta < self.convergence_tolerance:
-                convergence += 1
-                self.logger.info("Updated convergence level to: %u" %
-                                 convergence)
-                if convergence == self.convergence_level:
-                    self.logger.info("Reached target convergence level.")
-                    break
-            else:
-                convergence = max(0, convergence-1)
-
-            if self.iteration_limit is not None:
-                if iteration_number == self.iteration_limit:
-                    self.logger.warn("Reached iteration limit. Stopping.")
-                    break
-
-            iteration_number += 1
-
-        return energy, convergence
+            status = self._controller.check(energy)