add reference L_BFGS implementation

nifty4/minimization/__init__.py

@@ -8,14 +8,16 @@ from .nonlinear_cg import NonlinearCG
 from .descent_minimizer import DescentMinimizer
 from .steepest_descent import SteepestDescent
 from .vl_bfgs import VL_BFGS
+from .l_bfgs import L_BFGS
 from .relaxed_newton import RelaxedNewton
 from .scipy_minimizer import ScipyMinimizer, NewtonCG, L_BFGS_B
 from .energy import Energy
 from .quadratic_energy import QuadraticEnergy
 from .line_energy import LineEnergy
 
-__all__ = ["LineSearch", "LineSearchStrongWolfe", "IterationController",
-           "GradientNormController",
+__all__ = ["LineSearch", "LineSearchStrongWolfe",
+           "IterationController", "GradientNormController",
            "Minimizer", "ConjugateGradient", "NonlinearCG", "DescentMinimizer",
            "SteepestDescent", "VL_BFGS", "RelaxedNewton", "ScipyMinimizer",
-           "NewtonCG", "L_BFGS_B", "Energy", "QuadraticEnergy", "LineEnergy"]
+           "NewtonCG", "L_BFGS_B", "Energy", "QuadraticEnergy", "LineEnergy",
+           "L_BFGS"]

nifty4/minimization/l_bfgs.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-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 builtins import range
+import numpy as np
+from .descent_minimizer import DescentMinimizer
+from .line_search_strong_wolfe import LineSearchStrongWolfe
+from .. import dobj
+
+
+class L_BFGS(DescentMinimizer):
+
+    def __init__(self, controller, line_searcher=LineSearchStrongWolfe(),
+                 max_history_length=5):
+        super(L_BFGS, self).__init__(controller=controller,
+                                     line_searcher=line_searcher)
+        self.max_history_length = max_history_length
+
+    def __call__(self, energy):
+        self.reset()
+        return super(L_BFGS, self).__call__(energy)
+
+    def reset(self):
+        self._k = 0
+        self._s = [None]*self.max_history_length
+        self._y = [None]*self.max_history_length
+
+    def get_descent_direction(self, energy):
+        x = energy.position
+        s = self._s
+        y = self._y
+        k = self._k
+        maxhist = self.max_history_length
+        gradient = energy.gradient
+
+        nhist = min(k, maxhist)
+        alpha = [None]*maxhist
+        p = -gradient
+        if k > 0:
+            idx = (k-1) % maxhist
+            s[idx] = x-self._lastx
+            y[idx] = gradient-self._lastgrad
+        if nhist > 0:
+            for i in range(k-1, k-nhist-1, -1):
+                idx = i % maxhist
+                alpha[idx] = s[idx].vdot(p)/s[idx].vdot(y[idx])
+                p -= alpha[idx]*y[idx]
+            idx = (k-1) % maxhist
+            fact = s[idx].vdot(y[idx]) / y[idx].vdot(y[idx])
+            if fact <= 0.:
+                dobj.mprint("L-BFGS curvature not positive definite!")
+            p *= fact
+            for i in range(k-nhist, k):
+                idx = i % maxhist
+                beta = y[idx].vdot(p) / s[idx].vdot(y[idx])
+                p += (alpha[idx]-beta)*s[idx]
+        self._lastx = x
+        self._lastgrad = gradient
+        self._k += 1
+        return p

test/test_minimization/test_minimizers.py

@@ -34,7 +34,8 @@ minimizers = ['ift.VL_BFGS(IC)',
               'ift.NonlinearCG(IC, "Fletcher-Reeves")',
               'ift.NonlinearCG(IC, "5.49")',
               'ift.NewtonCG(IC)',
-              'ift.L_BFGS_B(IC)']
+              'ift.L_BFGS_B(IC)',
+              'ift.L_BFGS(IC)']
 
 newton_minimizers = ['ift.RelaxedNewton(IC)']
 quadratic_only_minimizers = ['ift.ConjugateGradient(IC)']