Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

Merge pull request #127 from dlaxalde/enh/optimize/ncg-hess

ENH: optimize/minimize: combine hess and hessp parameters
  • Loading branch information...
commit c9c2d66701583984589c88c08c478e2fc6d9f4ec 2 parents e5185f1 + c5f12ac
@teoliphant teoliphant authored
View
20 scipy/optimize/minimize.py
@@ -19,7 +19,7 @@
def minimize(fun, x0, args=(), method='Nelder-Mead', jac=None, hess=None,
- hessp=None, options=dict(), full_output=False, callback=None,
+ options=dict(), full_output=False, callback=None,
retall=False):
"""
Minimization of scalar function of one or more variables.
@@ -39,15 +39,13 @@ def minimize(fun, x0, args=(), method='Nelder-Mead', jac=None, hess=None,
jac : callable, optional
Jacobian of objective function (if None, Jacobian will be
estimated numerically). Only for CG, BFGS, Newton-CG.
- hess, hessp : callable, optional
- Hessian of objective function or Hessian of objective function
- times an arbitrary vector p. Only for Newton-CG.
- Only one of `hessp` or `hess` needs to be given. If `hess` is
- provided, then `hessp` will be ignored. If neither `hess` nor
- `hessp` is provided, then the hessian product will be approximated
- using finite differences on `jac`. `hessp` must compute the hessian
- times an arbitrary vector. If it is not given, finite-differences
- on `jac` are used to compute it.
+ hess : callable, optional
+ Hessian of objective function. Only for Newton-CG.
+ The function `hess` can either return the Hessian matrix of `fun`
+ or the Hessian matrix times an arbitrary vector, in which case
+ it accepts an extra argument `p` as `hess(x, p, *args)`.
+ If `hess` is None, the Hessian will be approximated using finite
+ differences on `jac`.
options : dict, optional
A dictionary of solver options. All methods accept the following
generic options:
@@ -203,7 +201,7 @@ def minimize(fun, x0, args=(), method='Nelder-Mead', jac=None, hess=None,
return _minimize_bfgs(fun, x0, args, jac, options, full_output,
retall, callback)
elif method.lower() == 'newton-cg':
- return _minimize_newtoncg(fun, x0, args, jac, hess, hessp, options,
+ return _minimize_newtoncg(fun, x0, args, jac, hess, options,
full_output, retall, callback)
elif method.lower() == 'anneal':
if callback:
View
33 scipy/optimize/optimize.py
@@ -28,6 +28,7 @@
from linesearch import \
line_search_BFGS, line_search_wolfe1, line_search_wolfe2, \
line_search_wolfe2 as line_search
+import inspect
# standard status messages of optimizers
@@ -1148,13 +1149,20 @@ def fmin_ncg(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5,
'maxiter': maxiter,
'disp': disp}
+ if fhess is not None:
+ hess = fhess
+ elif fhess_p is not None:
+ hess = fhess_p
+ else:
+ hess = None
+
# force full_output if retall=True to preserve backwards compatibility
if retall and not full_output:
- out = _minimize_newtoncg(f, x0, args, fprime, fhess, fhess_p, opts,
+ out = _minimize_newtoncg(f, x0, args, fprime, hess, opts,
full_output=True, retall=retall,
callback=callback)
else:
- out = _minimize_newtoncg(f, x0, args, fprime, fhess, fhess_p, opts,
+ out = _minimize_newtoncg(f, x0, args, fprime, hess, opts,
full_output, retall, callback)
if full_output:
@@ -1171,7 +1179,7 @@ def fmin_ncg(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5,
else:
return out
-def _minimize_newtoncg(fun, x0, args=(), jac=None, hess=None, hessp=None,
+def _minimize_newtoncg(fun, x0, args=(), jac=None, hess=None,
options={}, full_output=0, retall=0, callback=None):
"""
Minimization of scalar function of one or more variables using the
@@ -1197,8 +1205,23 @@ def _minimize_newtoncg(fun, x0, args=(), jac=None, hess=None, hessp=None,
raise ValueError('Jacobian is required for Newton-CG method')
f = fun
fprime = jac
- fhess_p = hessp
- fhess = hess
+ if hess is None:
+ fhess = None
+ fhess_p = None
+ else:
+ # check hessian type based on the number of arguments
+ fun_args = inspect.getargspec(fun)[0]
+ hess_args = inspect.getargspec(hess)[0]
+ if len(hess_args) == len(fun_args):
+ fhess = hess
+ fhess_p = None
+ elif len(hess_args) == len(fun_args) + 1:
+ fhess = None
+ fhess_p = hess
+ else:
+ raise ValueError('The number of arguments of the Hessian '
+ 'function does not agree with that of the '
+ 'objective function.')
# retrieve useful options
avextol = options.get('xtol', 1e-5)
epsilon = options.get('eps', _epsilon)
View
128 scipy/optimize/tests/test_optimize.py
@@ -14,10 +14,8 @@
assert_equal, assert_, TestCase, run_module_suite
from scipy import optimize
-from numpy import array, zeros, float64, dot, log, exp, inf, sin, cos
import numpy as np
from scipy.optimize.tnc import RCSTRINGS, MSG_NONE
-import numpy.random
from math import pow
class TestOptimize(TestCase):
@@ -26,10 +24,10 @@ class TestOptimize(TestCase):
Computational Linguistics, vol 22, num 1, pp 39--72, 1996.)
"""
def setUp(self):
- self.F = array([[1,1,1],[1,1,0],[1,0,1],[1,0,0],[1,0,0]])
- self.K = array([1., 0.3, 0.5])
- self.startparams = zeros(3, float64)
- self.solution = array([0., -0.524869316, 0.487525860])
+ self.F = np.array([[1,1,1],[1,1,0],[1,0,1],[1,0,0],[1,0,0]])
+ self.K = np.array([1., 0.3, 0.5])
+ self.startparams = np.zeros(3, np.float64)
+ self.solution = np.array([0., -0.524869316, 0.487525860])
self.maxiter = 1000
self.funccalls = 0
self.gradcalls = 0
@@ -40,19 +38,31 @@ def func(self, x):
self.funccalls += 1
if self.funccalls > 6000:
raise RuntimeError("too many iterations in optimization routine")
- log_pdot = dot(self.F, x)
- logZ = log(sum(exp(log_pdot)))
- f = logZ - dot(self.K, x)
+ log_pdot = np.dot(self.F, x)
+ logZ = np.log(sum(np.exp(log_pdot)))
+ f = logZ - np.dot(self.K, x)
self.trace.append(x)
return f
def grad(self, x):
self.gradcalls += 1
- log_pdot = dot(self.F, x)
- logZ = log(sum(exp(log_pdot)))
- p = exp(log_pdot - logZ)
- return dot(self.F.transpose(), p) - self.K
+ log_pdot = np.dot(self.F, x)
+ logZ = np.log(sum(np.exp(log_pdot)))
+ p = np.exp(log_pdot - logZ)
+ return np.dot(self.F.transpose(), p) - self.K
+
+
+ def hess(self, x):
+ log_pdot = np.dot(self.F, x)
+ logZ = np.log(sum(np.exp(log_pdot)))
+ p = np.exp(log_pdot - logZ)
+ return np.dot(self.F.T,
+ np.dot(np.diag(p), self.F - np.dot(self.F.T, p)))
+
+
+ def hessp(self, x, p):
+ return np.dot(self.hess(x), p)
def test_cg(self, use_wrapper=False):
@@ -280,6 +290,76 @@ def test_ncg(self, use_wrapper=False):
[-4.35700753e-07, -5.24869401e-01, 4.87527774e-01]],
atol=1e-6, rtol=1e-7), self.trace[:5])
+ def test_ncg_hess(self, use_wrapper=False):
+ """ Newton conjugate gradient with Hessian """
+ if use_wrapper:
+ opts = {'maxit': self.maxiter, 'disp': False}
+ retval = optimize.minimize(self.func, self.startparams,
+ method='Newton-CG', jac=self.grad,
+ hess = self.hess,
+ args=(), options=opts,
+ full_output=False, retall=False)
+ else:
+ retval = optimize.fmin_ncg(self.func, self.startparams, self.grad,
+ fhess = self.hess,
+ args=(), maxiter=self.maxiter,
+ full_output=False, disp=False,
+ retall=False)
+
+ params = retval
+
+ err = abs(self.func(params) - self.func(self.solution))
+ #print "NCG: Difference is: " + str(err)
+ assert_(err < 1e-6)
+
+ # Ensure that function call counts are 'known good'; these are from
+ # Scipy 0.7.0. Don't allow them to increase.
+ assert_(self.funccalls == 7, self.funccalls)
+ assert_(self.gradcalls <= 18, self.gradcalls) # 0.9.0
+ #assert_(self.gradcalls == 18, self.gradcalls) # 0.8.0
+ #assert_(self.gradcalls == 22, self.gradcalls) # 0.7.0
+
+ # Ensure that the function behaves the same; this is from Scipy 0.7.0
+ assert_(np.allclose(self.trace[3:5],
+ [[-4.35700753e-07, -5.24869435e-01, 4.87527480e-01],
+ [-4.35700753e-07, -5.24869401e-01, 4.87527774e-01]],
+ atol=1e-6, rtol=1e-7), self.trace[:5])
+
+ def test_ncg_hessp(self, use_wrapper=False):
+ """ Newton conjugate gradient with Hessian times a vector p """
+ if use_wrapper:
+ opts = {'maxit': self.maxiter, 'disp': False}
+ retval = optimize.minimize(self.func, self.startparams,
+ method='Newton-CG', jac=self.grad,
+ hess = self.hessp,
+ args=(), options=opts,
+ full_output=False, retall=False)
+ else:
+ retval = optimize.fmin_ncg(self.func, self.startparams, self.grad,
+ fhess_p = self.hessp,
+ args=(), maxiter=self.maxiter,
+ full_output=False, disp=False,
+ retall=False)
+
+ params = retval
+
+ err = abs(self.func(params) - self.func(self.solution))
+ #print "NCG: Difference is: " + str(err)
+ assert_(err < 1e-6)
+
+ # Ensure that function call counts are 'known good'; these are from
+ # Scipy 0.7.0. Don't allow them to increase.
+ assert_(self.funccalls == 7, self.funccalls)
+ assert_(self.gradcalls <= 18, self.gradcalls) # 0.9.0
+ #assert_(self.gradcalls == 18, self.gradcalls) # 0.8.0
+ #assert_(self.gradcalls == 22, self.gradcalls) # 0.7.0
+
+ # Ensure that the function behaves the same; this is from Scipy 0.7.0
+ assert_(np.allclose(self.trace[3:5],
+ [[-4.35700753e-07, -5.24869435e-01, 4.87527480e-01],
+ [-4.35700753e-07, -5.24869401e-01, 4.87527774e-01]],
+ atol=1e-6, rtol=1e-7), self.trace[:5])
+
def test_l_bfgs_b(self, use_wrapper=False):
""" limited-memory bound-constrained BFGS algorithm
@@ -324,6 +404,10 @@ def test_minimize(self):
self.setUp()
self.test_ncg(True)
self.setUp()
+ self.test_ncg_hess(True)
+ self.setUp()
+ self.test_ncg_hessp(True)
+ self.setUp()
self.test_neldermead(True)
self.setUp()
self.test_powell(True)
@@ -351,7 +435,7 @@ def test_fminbound(self):
x = optimize.fminbound(lambda x: (x - 1.5)**2 - 0.8, 1, 5)
assert_(abs(x - 1.5) < 1e-6)
x = optimize.fminbound(lambda x: (x - 1.5)**2 - 0.8,
- numpy.array([1]), numpy.array([5]))
+ np.array([1]), np.array([5]))
assert_(abs(x - 1.5) < 1e-6)
assert_raises(ValueError,
optimize.fminbound, lambda x: (x - 1.5)**2 - 0.8, 5, 1)
@@ -384,7 +468,7 @@ def test1fg(x):
dif[1] = 200.0*(x[1]-pow(x[0],2))
dif[0] = -2.0*(x[0]*(dif[1]-1.0)+1.0)
return f, dif
- self.tests.append((test1fg, [-2,1], ([-inf,None],[-1.5,None]),
+ self.tests.append((test1fg, [-2,1], ([-np.inf,None],[-1.5,None]),
[1,1]))
def test2fg(x):
f = 100.0*pow((x[1]-pow(x[0],2)),2)+pow(1.0-x[0],2)
@@ -392,7 +476,7 @@ def test2fg(x):
dif[1] = 200.0*(x[1]-pow(x[0],2))
dif[0] = -2.0*(x[0]*(dif[1]-1.0)+1.0)
return f, dif
- self.tests.append((test2fg, [-2,1], [(-inf,None),(1.5,None)],
+ self.tests.append((test2fg, [-2,1], [(-np.inf,None),(1.5,None)],
[-1.2210262419616387,1.5]))
def test3fg(x):
@@ -401,7 +485,7 @@ def test3fg(x):
dif[0] = -2.0*(x[1]-x[0])*1.0e-5
dif[1] = 1.0-dif[0]
return f, dif
- self.tests.append((test3fg, [10,1], [(-inf,None),(0.0, None)],
+ self.tests.append((test3fg, [10,1], [(-np.inf,None),(0.0, None)],
[0,0]))
def test4fg(x):
@@ -414,9 +498,9 @@ def test4fg(x):
[1,0]))
def test5fg(x):
- f = sin(x[0]+x[1])+pow(x[0]-x[1],2)-1.5*x[0]+2.5*x[1]+1.0
+ f = np.sin(x[0]+x[1])+pow(x[0]-x[1],2)-1.5*x[0]+2.5*x[1]+1.0
dif = [0,0]
- v1 = cos(x[0]+x[1])
+ v1 = np.cos(x[0]+x[1])
v2 = 2.0*(x[0]-x[1])
dif[0] = v1+v2-1.5
@@ -439,7 +523,7 @@ def test38fg(x):
dif[3] = (180.0*(x[3]-pow(x[2],2))+20.2\
*(x[3]-1.0)+19.8*(x[1]-1.0))*1.0e-5
return f, dif
- self.tests.append((test38fg, array([-3,-1,-3,-1]), [(-10,10)]*4, [1]*4))
+ self.tests.append((test38fg, np.array([-3,-1,-3,-1]), [(-10,10)]*4, [1]*4))
def test45fg(x):
f = 2.0-x[0]*x[1]*x[2]*x[3]*x[4]/120.0
@@ -470,8 +554,8 @@ class TestRosen(TestCase):
def test_hess(self):
"""Compare rosen_hess(x) times p with rosen_hess_prod(x,p) (ticket #1248)"""
- x = array([3, 4, 5])
- p = array([2, 2, 2])
+ x = np.array([3, 4, 5])
+ p = np.array([2, 2, 2])
hp = optimize.rosen_hess_prod(x, p)
dothp = np.dot(optimize.rosen_hess(x), p)
assert_equal(hp, dothp)
Please sign in to comment.
Something went wrong with that request. Please try again.