Merge branch 'improve_minimization' into 'NIFTy_4'

Improve minimization

See merge request ift/NIFTy!230

nifty4/minimization/descent_minimizer.py

@@ -81,11 +81,13 @@ class DescentMinimizer(Minimizer):
 
             # compute a step length that reduces energy.value sufficiently
             try:
-                new_energy = self.line_searcher.perform_line_search(
+                new_energy, success = self.line_searcher.perform_line_search(
                     energy=energy, pk=self.get_descent_direction(energy),
                     f_k_minus_1=f_k_minus_1)
             except ValueError:
                 return energy, controller.ERROR
+            if not success:
+                self.reset()
 
             f_k_minus_1 = energy.value
 
@@ -103,6 +105,9 @@ class DescentMinimizer(Minimizer):
             if status != controller.CONTINUE:
                 return energy, status
 
+    def reset(self):
+        pass
+
     @abc.abstractmethod
     def get_descent_direction(self, energy):
         """ Calculates the next descent direction.

nifty4/minimization/energy.py

@@ -109,3 +109,19 @@ class Energy(NiftyMetaBase()):
             curvature of the potential at the given `position`.
         """
         raise NotImplementedError
+
+    def longest_step(self, dir):
+        """Returns the longest allowed step size along `dir`
+
+        Parameters
+        ----------
+        dir : Field
+            the search direction
+
+        Returns
+        -------
+        float or None
+            the longest allowed step when starting from `self.position` along
+            `dir`. If None, the step size is not limited.
+        """
+        return None

nifty4/minimization/line_search.py

@@ -54,5 +54,7 @@ class LineSearch(NiftyMetaBase()):
         -------
         Energy
             The new Energy object on the new position.
+        bool
+            whether the line search was considered successful or not
         """
         raise NotImplementedError

nifty4/minimization/line_search_strong_wolfe.py

@@ -41,7 +41,7 @@ class LineSearchStrongWolfe(LineSearch):
         Parameter for curvature condition rule. (Default: 0.9)
     max_step_size : float
         Maximum step allowed in to be made in the descent direction.
-        (Default: 1000000000)
+        (Default: 1e30)
     max_iterations : int, optional
         Maximum number of iterations performed by the line search algorithm.
         (Default: 100)
@@ -51,7 +51,7 @@ class LineSearchStrongWolfe(LineSearch):
     """
 
     def __init__(self, preferred_initial_step_size=None, c1=1e-4, c2=0.9,
-                 max_step_size=1000000000, max_iterations=100,
+                 max_step_size=1e30, max_iterations=100,
                  max_zoom_iterations=100):
 
         super(LineSearchStrongWolfe, self).__init__(
@@ -85,19 +85,26 @@ class LineSearchStrongWolfe(LineSearch):
         -------
         Energy
             The new Energy object on the new position.
+        bool
+            whether the line search was considered successful or not
         """
         le_0 = LineEnergy(0., energy, pk, 0.)
 
+        maxstepsize = energy.longest_step(pk)
+        if maxstepsize is None:
+            maxstepsize = self.max_step_size
+        maxstepsize = min(maxstepsize, self.max_step_size)
+
         # initialize the zero phis
         old_phi_0 = f_k_minus_1
         phi_0 = le_0.value
         phiprime_0 = le_0.directional_derivative
         if phiprime_0 == 0:
             dobj.mprint("Directional derivative is zero; assuming convergence")
-            return energy
+            return energy, False
         if phiprime_0 > 0:
             dobj.mprint("Error: search direction is not a descent direction")
-            raise ValueError("search direction must be a descent direction")
+            return energy, False
 
         # set alphas
         alpha0 = 0.
@@ -112,49 +119,44 @@ class LineSearchStrongWolfe(LineSearch):
                 alpha1 = 1.0
         else:
             alpha1 = 1.0/pk.norm()
+        alpha1 = min(alpha1, 0.99*maxstepsize)
 
         # start the minimization loop
         iteration_number = 0
         while iteration_number < self.max_iterations:
             iteration_number += 1
             if alpha1 == 0:
-                result_energy = le_0.energy
-                break
+                return le_0.energy, False
 
             le_alpha1 = le_0.at(alpha1)
             phi_alpha1 = le_alpha1.value
 
             if (phi_alpha1 > phi_0 + self.c1*alpha1*phiprime_0) or \
                ((phi_alpha1 >= phi_alpha0) and (iteration_number > 1)):
-                le_star = self._zoom(alpha0, alpha1, phi_0, phiprime_0,
-                                     phi_alpha0, phiprime_alpha0, phi_alpha1,
-                                     le_0)
-                result_energy = le_star.energy
-                break
+                return self._zoom(alpha0, alpha1, phi_0, phiprime_0,
+                                  phi_alpha0, phiprime_alpha0, phi_alpha1,
+                                  le_0)
 
             phiprime_alpha1 = le_alpha1.directional_derivative
             if abs(phiprime_alpha1) <= -self.c2*phiprime_0:
-                result_energy = le_alpha1.energy
-                break
+                return le_alpha1.energy, True
 
             if phiprime_alpha1 >= 0:
-                le_star = self._zoom(alpha1, alpha0, phi_0, phiprime_0,
-                                     phi_alpha1, phiprime_alpha1, phi_alpha0,
-                                     le_0)
-                result_energy = le_star.energy
-                break
+                return self._zoom(alpha1, alpha0, phi_0, phiprime_0,
+                                  phi_alpha1, phiprime_alpha1, phi_alpha0,
+                                  le_0)
 
             # update alphas
-            alpha0, alpha1 = alpha1, min(2*alpha1, self.max_step_size)
-            if alpha1 == self.max_step_size:
-                return le_alpha1.energy
+            alpha0, alpha1 = alpha1, min(2*alpha1, maxstepsize)
+            if alpha1 == maxstepsize:
+                dobj.mprint("max step size reached")
+                return le_alpha1.energy, False
 
             phi_alpha0 = phi_alpha1
             phiprime_alpha0 = phiprime_alpha1
-        else:
-            dobj.mprint("max iterations reached")
-            return le_alpha1.energy
-        return result_energy
+
+        dobj.mprint("max iterations reached")
+        return le_alpha1.energy, False
 
     def _zoom(self, alpha_lo, alpha_hi, phi_0, phiprime_0,
               phi_lo, phiprime_lo, phi_hi, le_0):
@@ -238,7 +240,7 @@ class LineSearchStrongWolfe(LineSearch):
                 phiprime_alphaj = le_alphaj.directional_derivative
                 # If the second Wolfe condition is met, return the result
                 if abs(phiprime_alphaj) <= -self.c2*phiprime_0:
-                    return le_alphaj
+                    return le_alphaj.energy, True
                 # If not, check the sign of the slope
                 if phiprime_alphaj*delta_alpha >= 0:
                     alpha_recent, phi_recent = alpha_hi, phi_hi
@@ -251,7 +253,7 @@ class LineSearchStrongWolfe(LineSearch):
 
         else:
             dobj.mprint("The line search algorithm (zoom) did not converge.")
-            return le_alphaj
+            return le_alphaj.energy, False
 
     def _cubicmin(self, a, fa, fpa, b, fb, c, fc):
         """Estimating the minimum with cubic interpolation.

nifty4/minimization/nonlinear_cg.py

@@ -56,8 +56,10 @@ class NonlinearCG(Minimizer):
         while True:
             grad_old = energy.gradient
             f_k = energy.value
-            energy = self._line_searcher.perform_line_search(energy, p,
-                                                             f_k_minus_1)
+            energy, success = self._line_searcher.perform_line_search(
+                energy, p, f_k_minus_1)
+            if not success:
+                return energy, controller.ERROR
             f_k_minus_1 = f_k
             status = self._controller.check(energy)
             if status != controller.CONTINUE:

nifty4/minimization/scipy_minimizer.py

@@ -16,7 +16,7 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
-from __future__ import division, print_function
+from __future__ import division
 from .minimizer import Minimizer
 from ..field import Field
 from .. import dobj
@@ -97,9 +97,9 @@ class ScipyMinimizer(Minimizer):
             status = self._controller.check(hlp._energy)
             return hlp._energy, self._controller.check(hlp._energy)
         if not r.success:
-            print("Problem in Scipy minimization:", r.message)
+            dobj.mprint("Problem in Scipy minimization:", r.message)
         else:
-            print("Problem in Scipy minimization")
+            dobj.mprint("Problem in Scipy minimization")
         return hlp._energy, self._controller.ERROR
 
 

nifty4/minimization/vl_bfgs.py

@@ -49,6 +49,9 @@ class VL_BFGS(DescentMinimizer):
         self._information_store = None
         return super(VL_BFGS, self).__call__(energy)
 
+    def reset(self):
+        self._information_store = None
+
     def get_descent_direction(self, energy):
         x = energy.position
         gradient = energy.gradient