Merge branch 'fix_nonlinear_cg' into 'NIFTy_4'

Improve nonlinear conjugate gradient minimizers

See merge request ift/NIFTy!231

nifty4/minimization/nonlinear_cg.py

@@ -25,14 +25,17 @@ class NonlinearCG(Minimizer):
     """ Nonlinear Conjugate Gradient scheme according to Polak-Ribiere.
 
     Algorithm 5.4 from Nocedal & Wright.
-    Eq. (5.41a) has been replaced by eq. (5.49)
 
     Parameters
     ----------
     controller : IterationController
         Object that decides when to terminate the minimization.
-    line_searcher : LineSearch, optional
-        The line search algorithm to be used
+    beta_heuristics : str
+        One of 'Polak-Ribiere', 'Fletcher-Reeves', 'Hestenes-Stiefel' or '5.49'
+
+    Notes
+    -----
+    No restarting procedure has been implemented yet.
 
     References
     ----------
@@ -40,9 +43,15 @@ class NonlinearCG(Minimizer):
     2006, Springer-Verlag New York
     """
 
-    def __init__(self, controller, line_searcher=LineSearchStrongWolfe()):
+    def __init__(self, controller, beta_heuristics='Polak-Ribiere'):
+        valid_beta_heuristics = ['Polak-Ribiere', 'Fletcher-Reeves',
+                                 'Hestenes-Stiefel', "5.49"]
+        if not (beta_heuristics in valid_beta_heuristics):
+            raise ValueError("beta heuristics must be either 'Polak-Ribiere', "
+                             "'Fletcher-Reeves', 'Hestenes-Stiefel, or '5.49'")
+        self._beta_heuristic = beta_heuristics
         self._controller = controller
-        self._line_searcher = line_searcher
+        self._line_searcher = LineSearchStrongWolfe(c2=0.1)
 
     def __call__(self, energy):
         controller = self._controller
@@ -65,6 +74,20 @@ class NonlinearCG(Minimizer):
             if status != controller.CONTINUE:
                 return energy, status
             grad_new = energy.gradient
-            gnnew = energy.gradient_norm
-            beta = gnnew*gnnew/(grad_new-grad_old).vdot(p).real
+
+            if self._beta_heuristic == 'Hestenes-Stiefel':
+                # Eq. (5.46) in Nocedal & Wright.
+                beta = max(0.0, (grad_new.vdot(grad_new-grad_old) /
+                                 (grad_new-grad_old).vdot(p)).real)
+            elif self._beta_heuristic == 'Polak-Ribiere':
+                # Eq. (5.44) in Nocedal & Wright. (with (5.45) additionally)
+                beta = max(0.0, (grad_new.vdot(grad_new-grad_old) /
+                                 (grad_old.vdot(grad_old))).real)
+            elif self._beta_heuristic == 'Fletcher-Reeves':
+                # Eq. (5.41a) in Nocedal & Wright.
+                beta = (grad_new.vdot(grad_new)/(grad_old.vdot(grad_old))).real
+            else:
+                # Eq. (5.49) in Nocedal & Wright.
+                beta = (grad_new.vdot(grad_new) /
+                        ((grad_new-grad_old).vdot(p))).real
             p = beta*p - grad_new

test/test_minimization/test_minimizers.py

@@ -24,22 +24,28 @@ from itertools import product
 from test.common import expand
 from nose.plugins.skip import SkipTest
 
+IC = ift.GradientNormController(tol_abs_gradnorm=1e-5, iteration_limit=1000)
+
 spaces = [ift.RGSpace([1024], distances=0.123), ift.HPSpace(32)]
-minimizers = [ift.SteepestDescent, ift.RelaxedNewton, ift.VL_BFGS,
-              ift.ConjugateGradient, ift.NonlinearCG,
-              ift.NewtonCG, ift.L_BFGS_B]
 
-minimizers2 = [ift.RelaxedNewton, ift.VL_BFGS, ift.NonlinearCG,
-               ift.NewtonCG, ift.L_BFGS_B]
+minimizers = ['ift.VL_BFGS(IC)',
+              'ift.NonlinearCG(IC, "Polak-Ribiere")',
+              # ift.NonlinearCG(IC, "Hestenes-Stiefel"),
+              'ift.NonlinearCG(IC, "Fletcher-Reeves")',
+              'ift.NonlinearCG(IC, "5.49")',
+              'ift.NewtonCG(IC)',
+              'ift.L_BFGS_B(IC)']
 
-minimizers3 = [ift.SteepestDescent, ift.RelaxedNewton, ift.VL_BFGS,
-               ift.NonlinearCG, ift.NewtonCG, ift.L_BFGS_B]
+newton_minimizers = ['ift.RelaxedNewton(IC)']
+quadratic_only_minimizers = ['ift.ConjugateGradient(IC)']
+slow_minimizers = ['ift.SteepestDescent(IC)']
 
 
 class Test_Minimizers(unittest.TestCase):
 
-    @expand(product(minimizers, spaces))
-    def test_quadratic_minimization(self, minimizer_class, space):
+    @expand(product(minimizers + newton_minimizers +
+                    quadratic_only_minimizers + slow_minimizers, spaces))
+    def test_quadratic_minimization(self, minimizer, space):
         np.random.seed(42)
         starting_point = ift.Field.from_random('normal', domain=space)*10
         covariance_diagonal = ift.Field.from_random(
@@ -47,10 +53,8 @@ class Test_Minimizers(unittest.TestCase):
         covariance = ift.DiagonalOperator(covariance_diagonal)
         required_result = ift.Field.ones(space, dtype=np.float64)
 
-        IC = ift.GradientNormController(tol_abs_gradnorm=1e-5,
-                                        iteration_limit=1000)
         try:
-            minimizer = minimizer_class(controller=IC)
+            minimizer = eval(minimizer)
             energy = ift.QuadraticEnergy(A=covariance, b=required_result,
                                          position=starting_point)
 
@@ -63,8 +67,8 @@ class Test_Minimizers(unittest.TestCase):
                         1./covariance_diagonal.to_global_data(),
                         rtol=1e-3, atol=1e-3)
 
-    @expand(product(minimizers2))
-    def test_rosenbrock(self, minimizer_class):
+    @expand(product(minimizers+newton_minimizers))
+    def test_rosenbrock(self, minimizer):
         try:
             from scipy.optimize import rosen, rosen_der, rosen_hess_prod
         except ImportError:
@@ -114,10 +118,8 @@ class Test_Minimizers(unittest.TestCase):
                 return ift.InversionEnabler(RBCurv(self._position),
                                             inverter=t2)
 
-        IC = ift.GradientNormController(tol_abs_gradnorm=1e-5,
-                                        iteration_limit=10000)
         try:
-            minimizer = minimizer_class(controller=IC)
+            minimizer = eval(minimizer)
             energy = RBEnergy(position=starting_point)
 
             (energy, convergence) = minimizer(energy)
@@ -128,10 +130,10 @@ class Test_Minimizers(unittest.TestCase):
         assert_allclose(energy.position.to_global_data(), 1.,
                         rtol=1e-3, atol=1e-3)
 
-    @expand(product(minimizers3))
-    def DISABLEDtest_nonlinear(self, minimizer_class):
+    @expand(product(minimizers+slow_minimizers))
+    def test_gauss(self, minimizer):
         space = ift.UnstructuredDomain((1,))
-        starting_point = ift.Field(space, val=5.)
+        starting_point = ift.Field(space, val=3.)
 
         class ExpEnergy(ift.Energy):
             def __init__(self, position):
@@ -154,10 +156,8 @@ class Test_Minimizers(unittest.TestCase):
                 return ift.DiagonalOperator(
                     ift.Field(self.position.domain, val=v))
 
-        IC = ift.GradientNormController(tol_abs_gradnorm=1e-10,
-                                        iteration_limit=10000)
         try:
-            minimizer = minimizer_class(controller=IC)
+            minimizer = eval(minimizer)
             energy = ExpEnergy(position=starting_point)
 
             (energy, convergence) = minimizer(energy)
@@ -167,3 +167,41 @@ class Test_Minimizers(unittest.TestCase):
         assert_equal(convergence, IC.CONVERGED)
         assert_allclose(energy.position.to_global_data(), 0.,
                         atol=1e-3)
+
+    @expand(product(minimizers+newton_minimizers+slow_minimizers))
+    def test_cosh(self, minimizer):
+        space = ift.UnstructuredDomain((1,))
+        starting_point = ift.Field(space, val=3.)
+
+        class CoshEnergy(ift.Energy):
+            def __init__(self, position):
+                super(CoshEnergy, self).__init__(position)
+
+            @property
+            def value(self):
+                x = self.position.to_global_data()[0]
+                return np.cosh(x)
+
+            @property
+            def gradient(self):
+                x = self.position.to_global_data()[0]
+                return ift.Field(self.position.domain, val=np.sinh(x))
+
+            @property
+            def curvature(self):
+                x = self.position.to_global_data()[0]
+                v = np.cosh(x)
+                return ift.DiagonalOperator(
+                    ift.Field(self.position.domain, val=v))
+
+        try:
+            minimizer = eval(minimizer)
+            energy = CoshEnergy(position=starting_point)
+
+            (energy, convergence) = minimizer(energy)
+        except NotImplementedError:
+            raise SkipTest
+
+        assert_equal(convergence, IC.CONVERGED)
+        assert_allclose(energy.position.to_global_data(), 0.,
+                        atol=1e-3)