implement a NIFTy-aware Newton-CG minimizer

demos/bernoulli_demo.py

@@ -66,15 +66,14 @@ if __name__ == '__main__':
     # Compute likelihood and Hamiltonian
     position = ift.from_random('normal', harmonic_space)
     likelihood = ift.BernoulliEnergy(p, data)
-    ic_cg = ift.GradientNormController(iteration_limit=50)
-    ic_newton = ift.GradientNormController(name='Newton', iteration_limit=30,
-                                           tol_abs_gradnorm=1e-3)
-    minimizer = ift.RelaxedNewton(ic_newton)
+    ic_newton = ift.DeltaEnergyController(name='Newton', iteration_limit=100,
+                                          tol_rel_deltaE=1e-8)
+    minimizer = ift.NewtonCG(ic_newton)
     ic_sampling = ift.GradientNormController(iteration_limit=100)
 
     # Minimize the Hamiltonian
     H = ift.Hamiltonian(likelihood, ic_sampling)
-    H = ift.EnergyAdapter(position, H, ic_cg)
+    H = ift.EnergyAdapter(position, H)
     # minimizer = ift.L_BFGS(ic_newton)
     H, convergence = minimizer(H)
 

demos/getting_started_2.py

@@ -87,14 +87,13 @@ if __name__ == '__main__':
 
     # Compute likelihood and Hamiltonian
     likelihood = ift.PoissonianEnergy(lamb, data)
-    ic_cg = ift.GradientNormController(iteration_limit=50)
-    ic_newton = ift.GradientNormController(name='Newton', iteration_limit=50,
-                                           tol_abs_gradnorm=1e-3)
-    minimizer = ift.RelaxedNewton(ic_newton)
+    ic_newton = ift.DeltaEnergyController(name='Newton', iteration_limit=100,
+                                          tol_rel_deltaE=1e-8)
+    minimizer = ift.NewtonCG(ic_newton)
 
     # Minimize the Hamiltonian
     H = ift.Hamiltonian(likelihood)
-    H = ift.EnergyAdapter(position, H, ic_cg)
+    H = ift.EnergyAdapter(position, H)
     H, convergence = minimizer(H)
 
     # Plot results

demos/getting_started_3.py

@@ -72,14 +72,10 @@ if __name__ == '__main__':
         mean=data, covariance=N)(signal_response)
 
     # set up minimization and inversion schemes
-    ic_cg = ift.GradientNormController(iteration_limit=10)
     ic_sampling = ift.GradientNormController(iteration_limit=100)
     ic_newton = ift.DeltaEnergyController(
         name='Newton', tol_rel_deltaE=1e-8, iteration_limit=100)
-    minimizer = ift.RelaxedNewton(ic_newton)
-    # minimizer = ift.VL_BFGS(ic_newton)
-    # minimizer = ift.NewtonCG(xtol=1e-10, maxiter=100, disp=True)
-    # minimizer = ift.L_BFGS_B(ftol=1e-10, gtol=1e-5, maxiter=100, maxcor=20, disp=True)
+    minimizer = ift.NewtonCG(ic_newton)
 
     # build model Hamiltonian
     H = ift.Hamiltonian(likelihood, ic_sampling)
@@ -100,7 +96,7 @@ if __name__ == '__main__':
                    for _ in range(N_samples)]
 
         KL = ift.SampledKullbachLeiblerDivergence(H, samples)
-        KL = ift.EnergyAdapter(position, KL, ic_cg, constants=["xi"])
+        KL = ift.EnergyAdapter(position, KL)
         KL, convergence = minimizer(KL)
         position = KL.position
 

demos/polynomial_fit.py

@@ -90,7 +90,7 @@ H = ift.Hamiltonian(likelihood, IC)
 H = ift.EnergyAdapter(params, H, IC)
 
 # Minimize
-minimizer = ift.RelaxedNewton(IC)
+minimizer = ift.NewtonCG(IC)
 H, _ = minimizer(H)
 
 # Draw posterior samples

nifty5/__init__.py

@@ -60,9 +60,9 @@ from .minimization.minimizer import Minimizer
 from .minimization.conjugate_gradient import ConjugateGradient
 from .minimization.nonlinear_cg import NonlinearCG
 from .minimization.descent_minimizers import (
-    DescentMinimizer, SteepestDescent, VL_BFGS, L_BFGS, RelaxedNewton)
-from .minimization.scipy_minimizer import (ScipyMinimizer, NewtonCG, L_BFGS_B,
-                                           ScipyCG)
+    DescentMinimizer, SteepestDescent, VL_BFGS, L_BFGS, RelaxedNewton,
+    NewtonCG)
+from .minimization.scipy_minimizer import (ScipyMinimizer, L_BFGS_B, ScipyCG)
 from .minimization.energy import Energy
 from .minimization.quadratic_energy import QuadraticEnergy
 from .minimization.line_energy import LineEnergy

nifty5/minimization/descent_minimizers.py

@@ -153,6 +153,57 @@ class RelaxedNewton(DescentMinimizer):
         return -energy.metric.inverse_times(energy.gradient)
 
 
+class NewtonCG(DescentMinimizer):
+    """ Calculates the descent direction according to a Newton-CG scheme.
+
+    Algorithm derived from SciPy sources.
+    """
+
+    def __init__(self, controller, line_searcher=None):
+        if line_searcher is None:
+            line_searcher = LineSearchStrongWolfe(
+                preferred_initial_step_size=1.)
+        super(NewtonCG, self).__init__(controller=controller,
+                                       line_searcher=line_searcher)
+
+    def get_descent_direction(self, energy):
+        float64eps = np.finfo(np.float64).eps
+        grad = energy.gradient
+        maggrad = abs(grad).sum()
+        eta = np.min([0.5, np.sqrt(maggrad)])
+        termcond = eta*maggrad
+        xsupi = energy.position*0
+        ri = grad
+        psupi = -ri
+        i = 0
+        dri0 = ri.vdot(ri)
+
+        while True:
+            if abs(ri).sum() <= termcond:
+                return xsupi
+            Ap = energy.metric(psupi)
+            # check curvature
+            curv = psupi.vdot(Ap)
+            if 0 <= curv <= 3 * float64eps:
+                return xsupi
+            elif curv < 0:
+                if (i > 0):
+                    return xsupi
+                else:
+                    return (dri0/curv) * grad  # fall back to steepest descent
+            alphai = dri0/curv
+            xsupi = xsupi + alphai*psupi
+            ri = ri + alphai*Ap
+            dri1 = ri.vdot(ri)
+            psupi = (dri1/dri0)*psupi - ri
+            i += 1
+            dri0 = dri1  # update numpy.dot(ri,ri) for next time.
+
+        # curvature keeps increasing, bail out
+        raise ValueError("Warning: CG iterations didn't converge. "
+                         "The Hessian is not positive definite.")
+
+
 class L_BFGS(DescentMinimizer):
     def __init__(self, controller, line_searcher=LineSearchStrongWolfe(),
                  max_history_length=5):

nifty5/minimization/scipy_minimizer.py

@@ -139,17 +139,6 @@ class ScipyMinimizer(Minimizer):
         return hlp._energy, IterationController.CONVERGED
 
 
-def NewtonCG(xtol, maxiter, disp=False):
-    """Returns a ScipyMinimizer object carrying out the Newton-CG algorithm.
-
-    See Also
-    --------
-    ScipyMinimizer
-    """
-    options = {"xtol": xtol, "maxiter": maxiter, "disp": disp}
-    return ScipyMinimizer("Newton-CG", options, True, None)
-
-
 def L_BFGS_B(ftol, gtol, maxiter, maxcor=10, disp=False, bounds=None):
     """Returns a ScipyMinimizer object carrying out the L-BFGS-B algorithm.
 

test/test_minimization/test_minimizers.py

@@ -34,9 +34,9 @@ minimizers = ['ift.VL_BFGS(IC)',
               # 'ift.NonlinearCG(IC, "Hestenes-Stiefel"),
               'ift.NonlinearCG(IC, "Fletcher-Reeves")',
               'ift.NonlinearCG(IC, "5.49")',
-              'ift.NewtonCG(xtol=1e-5, maxiter=1000)',
               'ift.L_BFGS_B(ftol=1e-10, gtol=1e-5, maxiter=1000)',
-              'ift.L_BFGS(IC)']
+              'ift.L_BFGS(IC)',
+              'ift.NewtonCG(IC)']
 
 newton_minimizers = ['ift.RelaxedNewton(IC)']
 quadratic_only_minimizers = ['ift.ConjugateGradient(IC)',