From 1b82593bcd883e4070928b4f0564a26eabda51a4 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Tue, 13 Mar 2018 13:37:30 +0100
Subject: [PATCH] add reference L_BFGS implementation
---
nifty4/minimization/__init__.py | 8 ++-
nifty4/minimization/l_bfgs.py | 76 +++++++++++++++++++++++
test/test_minimization/test_minimizers.py | 3 +-
3 files changed, 83 insertions(+), 4 deletions(-)
create mode 100644 nifty4/minimization/l_bfgs.py
diff --git a/nifty4/minimization/__init__.py b/nifty4/minimization/__init__.py
index 0dd11d63c..b24aa2768 100644
--- a/nifty4/minimization/__init__.py
+++ b/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"]
diff --git a/nifty4/minimization/l_bfgs.py b/nifty4/minimization/l_bfgs.py
new file mode 100644
index 000000000..c7e910e34
--- /dev/null
+++ b/nifty4/minimization/l_bfgs.py
@@ -0,0 +1,76 @@
+# 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
diff --git a/test/test_minimization/test_minimizers.py b/test/test_minimization/test_minimizers.py
index 76807ab92..222b91c7a 100644
--- a/test/test_minimization/test_minimizers.py
+++ b/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)']
--
GitLab