...
 
Commits (29)
# 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 print_function
import numpy as np
import nifty4 as ift
IC = ift.GradientNormController(
tol_abs_gradnorm=1e-6, iteration_limit=1000)
calls = 0
def mprint(*args):
if ift.dobj.rank == 0:
print(*args)
def test_rosenbrock_convex(minimizer):
np.random.seed(42)
space = ift.UnstructuredDomain((2,))
starting_point = ift.Field.from_random('normal', domain=space)*10
class RBEnergy(ift.Energy):
def __init__(self, position, a=1., b=100.):
super(RBEnergy, self).__init__(position)
self.a = a
self.b = b
global calls
calls += 1
@property
def value(self):
x, y = self.position.to_global_data()
return (self.a-x)*(self.a-x)+self.b*(y-x*x)*(y-x*x)
@property
def gradient(self):
x, y = self.position.to_global_data()
v0 = -2*(self.a-x)-4*self.b*x*(y-x*x)
v1 = 2*self.b*(y-x*x)
return ift.Field.from_global_data(space, np.array([v0, v1]))
@property
def curvature(self):
class RBCurv(ift.EndomorphicOperator):
def __init__(self, loc, a, b):
self._x, self._y = loc.to_global_data()
self.a = a
self.b = b
@property
def domain(self):
return space
@property
def capability(self):
return self.TIMES
def apply(self, x, mode):
x = x.to_global_data()
v0 = (2+self.b*8*self._x**2)*x[0]
v0 -= self.b*4*self._x*x[1]
v1 = -self.b*4*self._x*x[0]
v1 += 2*self.b*x[1]
global calls
calls += 1
return ift.Field.from_global_data(space,
np.array([v0, v1]))
t1 = ift.GradientNormController(tol_abs_gradnorm=1e-6,
iteration_limit=1000)
t2 = ift.ConjugateGradient(controller=t1)
return ift.InversionEnabler(RBCurv(self._position, self.a, self.b),
inverter=t2)
energy = RBEnergy(position=starting_point)
minimizer = minimizer(controller=IC)
energy, convergence = minimizer(energy)
return energy
def test_Ndim_rosenbrocklike_convex(minimizer, Ndim=3):
np.random.seed(42)
space = ift.UnstructuredDomain((Ndim,))
starting_point = ift.Field.from_random('normal', domain=space)*10
class RBLikeEnergy(ift.Energy):
def __init__(self, position, a=1., b=100.):
super(RBLikeEnergy, self).__init__(position)
self.a = a
self.b = b
global calls
calls += 1
@property
def value(self):
x = self.position.to_global_data()
t1 = self.a-x[0]
t2 = x[1:]-x[:-1]**3
return 0.5*t1**2+0.5*self.b*np.dot(t2, t2)
@property
def gradient(self):
res = np.zeros(space.shape)
x = self.position.to_global_data()
t1 = self.a-x[0]
t2 = x[1:]-x[:-1]**3
res[0] = -t1
res[1:] += self.b*t2
res[:-1] += -3*self.b*x[:-1]**2*t2
return ift.Field.from_global_data(space, res)
@property
def curvature(self):
class RBCurv(ift.EndomorphicOperator):
def __init__(self, loc, a, b):
self._x = loc.to_global_data()
self.a = a
self.b = b
@property
def domain(self):
return space
@property
def capability(self):
return self.TIMES
def apply(self, z, mode):
y = z.to_global_data()
res = np.zeros(space.shape)
x = self._x
res[0] = y[0]
dt2 = y[1:] - 3*x[:-1]**2*y[:-1]
res[1:] += self.b*dt2
res[:-1] += -self.b*3*x[:-1]**2*dt2
global calls
calls += 1
return ift.Field.from_global_data(space, res)
t1 = ift.GradientNormController(tol_abs_gradnorm=1e-6,
iteration_limit=1000)
t2 = ift.ConjugateGradient(controller=t1)
return ift.InversionEnabler(RBCurv(self._position, self.a, self.b),
inverter=t2)
energy = RBLikeEnergy(position=starting_point)
minimizer = minimizer(controller=IC)
energy, convergence = minimizer(energy)
return energy
def verbose_test(func, minimizer):
mprint("Testing", minimizer)
global calls
calls = 0
E = func(minimizer)
mprint("Used", calls, "calls.")
mprint("Final energy:", E.value)
if __name__ == "__main__":
mprint("Standard Rosenbrock function\n")
verbose_test(test_rosenbrock_convex, ift.Yango)
verbose_test(test_rosenbrock_convex, ift.RelaxedNewton)
verbose_test(test_rosenbrock_convex, ift.NonlinearCG)
verbose_test(test_rosenbrock_convex, ift.L_BFGS)
mprint("\nHigher-dimensional Rosenbrock function\n")
verbose_test(test_Ndim_rosenbrocklike_convex, ift.Yango)
verbose_test(test_Ndim_rosenbrocklike_convex, ift.RelaxedNewton)
verbose_test(test_Ndim_rosenbrocklike_convex, ift.NonlinearCG)
verbose_test(test_Ndim_rosenbrocklike_convex, ift.L_BFGS)
......@@ -14,6 +14,7 @@ from .scipy_minimizer import ScipyMinimizer, NewtonCG, L_BFGS_B, ScipyCG
from .energy import Energy
from .quadratic_energy import QuadraticEnergy
from .line_energy import LineEnergy
from .yango import Yango
from .energy_sum import EnergySum
__all__ = ["LineSearch", "LineSearchStrongWolfe",
......@@ -21,4 +22,4 @@ __all__ = ["LineSearch", "LineSearchStrongWolfe",
"Minimizer", "ConjugateGradient", "NonlinearCG", "DescentMinimizer",
"SteepestDescent", "VL_BFGS", "RelaxedNewton", "ScipyMinimizer",
"NewtonCG", "L_BFGS_B", "ScipyCG", "Energy", "QuadraticEnergy",
"LineEnergy", "L_BFGS", "EnergySum"]
"LineEnergy", "L_BFGS", "EnergySum", "Yango"]
# 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 .minimizer import Minimizer
from ..logger import logger
from .line_search_strong_wolfe import LineSearchStrongWolfe
import numpy as np
_default_LS = LineSearchStrongWolfe(c2=0.1, preferred_initial_step_size=1.)
class Yango(Minimizer):
""" Nonlinear conjugate gradient using curvature
The YANGO (Yet Another Nonlinear conjugate Gradient Optimizer)
uses the curvature to make estimates about suitable descent
directions. It takes the step that lets it go directly to
the second order minimum in the subspace spanned by the last
descent direction and the new gradient.
Parameters
----------
controller : IterationController
Object that decides when to terminate the minimization.
Notes
-----
No restarting procedure has been implemented yet.
"""
def __init__(self, controller, line_searcher=_default_LS):
self._controller = controller
self._line_searcher = line_searcher
def __call__(self, energy):
controller = self._controller
status = controller.start(energy)
if status != controller.CONTINUE:
return energy, status
f_k_minus_1 = None
p = -energy.gradient
A_k = energy.curvature
Ap = A_k(p)
pAp = p.vdot(Ap)
pp = p.vdot(p)
energy, success = self._line_searcher.perform_line_search(
energy, (pp/pAp)*p, f_k_minus_1)
if not success:
return energy, controller.ERROR
A_k = energy.curvature
while True:
r = -energy.gradient
f_k = energy.value
Ar = A_k(r)
Ap = A_k(p)
rAr = r.vdot(Ar)
pAp = p.vdot(Ap)
pAr = p.vdot(Ar)
rAp = np.conj(pAr)
rp = r.vdot(p)
rr = r.vdot(r)
if rr == 0 or rAr == 0:
logger.warning(
"Warning: gradient norm 0, assuming convergence!")
return energy, controller.CONVERGED
det = pAp*rAr-np.abs(rAp*pAr)
if det < 0:
if rAr < 0:
logger.error(
"Error: negative curvature ({})".format(rAr))
return energy, controller.ERROR
# Try 1D Newton Step
energy, success = self._line_searcher.perform_line_search(
energy, (rr/rAr)*r, f_k_minus_1)
else:
a = (rAr*rp - rAp*rr)/det
b = (pAp*rr - pAr*rp)/det
energy, success = self._line_searcher.perform_line_search(
energy, p*a + r*b, f_k_minus_1)
p = r - p*(pAr/pAp)
if not success:
return energy, controller.ERROR
f_k_minus_1 = f_k
status = self._controller.check(energy)
if status != controller.CONTINUE:
return energy, status
A_k = energy.curvature
......@@ -37,7 +37,7 @@ minimizers = ['ift.VL_BFGS(IC)',
'ift.L_BFGS_B(ftol=1e-10, gtol=1e-5, maxiter=1000)',
'ift.L_BFGS(IC)']
newton_minimizers = ['ift.RelaxedNewton(IC)']
newton_minimizers = ['ift.RelaxedNewton(IC)', 'ift.Yango(IC)']
quadratic_only_minimizers = ['ift.ConjugateGradient(IC)',
'ift.ScipyCG(tol=1e-5, maxiter=300)']
slow_minimizers = ['ift.SteepestDescent(IC)']
......@@ -69,7 +69,7 @@ class Test_Minimizers(unittest.TestCase):
1./covariance_diagonal.to_global_data(),
rtol=1e-3, atol=1e-3)
@expand(product(minimizers+newton_minimizers))
@expand(product(minimizers))
def test_rosenbrock(self, minimizer):
try:
from scipy.optimize import rosen, rosen_der, rosen_hess_prod
......@@ -132,6 +132,73 @@ class Test_Minimizers(unittest.TestCase):
assert_allclose(energy.position.to_global_data(), 1.,
rtol=1e-3, atol=1e-3)
@expand(product(minimizers+newton_minimizers))
def test_rosenbrock_convex(self, minimizer):
np.random.seed(42)
space = ift.UnstructuredDomain((2,))
starting_point = ift.Field.from_random('normal', domain=space)*10
class RBEnergy(ift.Energy):
def __init__(self, position, a=1., b=100.):
super(RBEnergy, self).__init__(position)
self.a = a
self.b = b
@property
def value(self):
x, y = self.position.to_global_data()
return (self.a-x)*(self.a-x)+self.b*(y-x*x)*(y-x*x)
@property
def gradient(self):
x, y = self.position.to_global_data()
r0 = -2*(self.a-x)-4*self.b*x*(y-x*x)
r1 = 2*self.b*(y-x*x)
return ift.Field.from_global_data(space, np.array([r0, r1]))
@property
def curvature(self):
class RBCurv(ift.EndomorphicOperator):
def __init__(self, loc, a, b):
self._x, self._y = loc.to_global_data()
self.a = a
self.b = b
@property
def domain(self):
return space
@property
def capability(self):
return self.TIMES
def apply(self, x, mode):
x = x.to_global_data()
r0 = (2+self.b*8*self._x**2)*x[0]
r0 -= self.b*4*self._x*x[1]
r1 = -self.b*4*self._x*x[0]
r1 += 2*self.b*x[1]
return ift.Field.from_global_data(space,
np.array([r0, r1]))
t1 = ift.GradientNormController(tol_abs_gradnorm=1e-5,
iteration_limit=1000)
t2 = ift.ConjugateGradient(controller=t1)
return ift.InversionEnabler(
RBCurv(self._position, self.a, self.b), inverter=t2)
energy = RBEnergy(position=starting_point)
try:
minimizer = eval(minimizer)
energy = RBEnergy(position=starting_point)
(energy, convergence) = minimizer(energy)
except NotImplementedError:
raise SkipTest
assert_equal(convergence, IC.CONVERGED)
assert_allclose(energy.position.to_global_data(), 1.,
rtol=1e-3, atol=1e-3)
@expand(product(minimizers+slow_minimizers))
def test_gauss(self, minimizer):
space = ift.UnstructuredDomain((1,))
......