From 09016f9fc4b2c87734201894b9eeb24eebb34d2f Mon Sep 17 00:00:00 2001 From: QBatista Date: Wed, 7 Nov 2018 09:13:00 +0400 Subject: [PATCH] Minor fixes --- quantecon/optimize/__init__.py | 2 +- ...ultivar_maximization.py => nelder_mead.py} | 84 +++++++++++-------- quantecon/optimize/tests/__init__.py | 0 ...st_multivar_max.py => test_nelder_mead.py} | 67 +++++++-------- 4 files changed, 86 insertions(+), 67 deletions(-) rename quantecon/optimize/{multivar_maximization.py => nelder_mead.py} (80%) create mode 100644 quantecon/optimize/tests/__init__.py rename quantecon/optimize/tests/{test_multivar_max.py => test_nelder_mead.py} (81%) diff --git a/quantecon/optimize/__init__.py b/quantecon/optimize/__init__.py index fb24c413a..ea049aa03 100644 --- a/quantecon/optimize/__init__.py +++ b/quantecon/optimize/__init__.py @@ -3,5 +3,5 @@ """ from .scalar_maximization import brent_max -from .multivar_maximization import maximize, nelder_mead_algorithm +from .nelder_mead import nelder_mead from .root_finding import newton, newton_halley, newton_secant, bisect, brentq diff --git a/quantecon/optimize/multivar_maximization.py b/quantecon/optimize/nelder_mead.py similarity index 80% rename from quantecon/optimize/multivar_maximization.py rename to quantecon/optimize/nelder_mead.py index a8b29c9e6..d0f7592a5 100644 --- a/quantecon/optimize/multivar_maximization.py +++ b/quantecon/optimize/nelder_mead.py @@ -1,5 +1,6 @@ """ -Implements the Nelder-Mead algorithm for maximizing a multivariate function +Implements the Nelder-Mead algorithm for maximizing a function with one or more +variables. """ @@ -11,13 +12,15 @@ @njit -def maximize(fun, x0, bounds=np.array([[], []]).T, args=(), tol_f=1e-10, - tol_x=1e-10, max_iter=1000): +def nelder_mead(fun, x0, bounds=np.array([[], []]).T, args=(), tol_f=1e-10, + tol_x=1e-10, max_iter=1000): """ .. highlight:: none - Maximize a scalar function with multiple variables using the Nelder-Mead - method. + Maximize a scalar-valued function with one or more variables using the + Nelder-Mead method. + + This function is JIT-compiled in `nopython` mode using Numba. Parameters ---------- @@ -25,24 +28,26 @@ def maximize(fun, x0, bounds=np.array([[], []]).T, args=(), tol_f=1e-10, The objective function to be maximized. `fun(x, *args) -> float` where x is an 1-D array with shape (n,) and args is a tuple of the - fixed parameters needed to completely specify the function. + fixed parameters needed to completely specify the function. This + function must be JIT-compiled in `nopython` mode using Numba. x0 : ndarray(float, ndim=1) Initial guess. Array of real elements of size (n,), where ‘n’ is the number of independent variables. bounds: ndarray(float, ndim=2), optional - Bounds for each variable for proposed solution. Sequence of (min, max) - pairs for each element in x. The default is used to specify no bounds. + Bounds for each variable for proposed solution, encoded as a sequence + of (min, max) pairs for each element in x. The default option is used + to specify no bounds on x. args : tuple, optional Extra arguments passed to the objective function. tol_f : scalar(float), optional(default=1e-10) - Tolerance for the range convergence test. + Tolerance to be used for the function value convergence test. tol_x : scalar(float), optional(default=1e-10) - Tolerance for the domain convergence test. + Tolerance to be used for the function domain convergence test. max_iter : scalar(float), optional(default=1000) The maximum number of allowed iterations. @@ -53,11 +58,11 @@ def maximize(fun, x0, bounds=np.array([[], []]).T, args=(), tol_f=1e-10, A namedtuple containing the following items: :: - "x" : Approximate solution - "fun" : Approximate local maximum - "success" : 1 if successfully terminated, 0 otherwise + "x" : Approximate local maximizer + "fun" : Approximate local maximum value + "success" : 1 if the algorithm successfully terminated, 0 otherwise "nit" : Number of iterations - "final_simplex" : The vertices of the final simplex + "final_simplex" : Vertices of the final simplex Examples -------- @@ -70,6 +75,16 @@ def maximize(fun, x0, bounds=np.array([[], []]).T, args=(), tol_f=1e-10, results(x=array([0.99999814, 0.99999756]), fun=1.6936258239463265e-10, success=True, nit=110) + Notes + -------- + This algorithm has a long history of successful use in applications, but it + will usually be slower than an algorithm that uses first or second + derivative information. In practice, it can have poor performance in + high-dimensional problems and is not robust to minimizing complicated + functions. Additionally, there currently is no complete theory describing + when the algorithm will successfully converge to the minimum, or how fast + it will if it does. + References ---------- @@ -93,38 +108,37 @@ def maximize(fun, x0, bounds=np.array([[], []]).T, args=(), tol_f=1e-10, .. [7] Chase Coleman's tutorial on Nelder Mead - .. [8] https://github.com/scipy/scipy/blob/v1.1.0/scipy/optimize/optimize.py + .. [8] SciPy's Nelder-Mead implementation """ - n = x0.size - vertices = _initialize_simplex(x0, bounds) - results = nelder_mead_algorithm(fun, vertices, bounds, args=args, - tol_f=tol_f, tol_x=tol_x, - max_iter=max_iter) + results = _nelder_mead_algorithm(fun, vertices, bounds, args=args, + tol_f=tol_f, tol_x=tol_x, + max_iter=max_iter) return results @njit -def nelder_mead_algorithm(fun, vertices, bounds=np.array([[], []]).T, - args=(), ρ=1., χ=2., γ=0.5, σ=0.5, tol_f=1e-8, - tol_x=1e-8, max_iter=1000): - +def _nelder_mead_algorithm(fun, vertices, bounds=np.array([[], []]).T, + args=(), ρ=1., χ=2., γ=0.5, σ=0.5, tol_f=1e-8, + tol_x=1e-8, max_iter=1000): """ .. highlight:: none Implements the Nelder-Mead algorithm described in Lagarias et al. (1998) - modified to maximize instead of minimizing. + modified to maximize instead of minimizing. JIT-compiled in `nopython` + mode using Numba. Parameters ---------- fun : callable - The objective function to be minimized. + The objective function to be maximized. `fun(x, *args) -> float` where x is an 1-D array with shape (n,) and args is a tuple of the - fixed parameters needed to completely specify the function. + fixed parameters needed to completely specify the function. This + function must be JIT-compiled in `nopython` mode using Numba. vertices : ndarray(float, ndim=2) Initial simplex with shape (n+1, n) to be modified in-place. @@ -145,10 +159,10 @@ def nelder_mead_algorithm(fun, vertices, bounds=np.array([[], []]).T, Shrinkage parameter. Must be strictly between 0 and 1. tol_f : scalar(float), optional(default=1e-10) - Tolerance for the range convergence test. + Tolerance to be used for the function value convergence test. tol_x : scalar(float), optional(default=1e-10) - Tolerance for the domain convergence test. + Tolerance to be used for the function domain convergence test. max_iter : scalar(float), optional(default=1000) The maximum number of allowed iterations. @@ -277,7 +291,8 @@ def nelder_mead_algorithm(fun, vertices, bounds=np.array([[], []]).T, @njit def _initialize_simplex(x0, bounds): """ - Generates an initial simplex for the Nelder-Mead method. + Generates an initial simplex for the Nelder-Mead method. JIT-compiled in + `nopython` mode using Numba. Parameters ---------- @@ -318,6 +333,7 @@ def _initialize_simplex(x0, bounds): def _check_params(ρ, χ, γ, σ, bounds, n): """ Checks whether the parameters for the Nelder-Mead algorithm are valid. + JIT-compiled in `nopython` mode using Numba. Parameters ---------- @@ -360,7 +376,8 @@ def _check_params(ρ, χ, γ, σ, bounds, n): @njit def _check_bounds(x, bounds): """ - Checks whether `x` is within `bounds`. + Checks whether `x` is within `bounds`. JIT-compiled in `nopython` mode + using Numba. Parameters ---------- @@ -387,7 +404,7 @@ def _check_bounds(x, bounds): def _neg_bounded_fun(fun, bounds, x, args=()): """ Wrapper for bounding and taking the negative of `fun` for the - Nelder-Mead algorithm. + Nelder-Mead algorithm. JIT-compiled in `nopython` mode using Numba. Parameters ---------- @@ -395,7 +412,8 @@ def _neg_bounded_fun(fun, bounds, x, args=()): The objective function to be minimized. `fun(x, *args) -> float` where x is an 1-D array with shape (n,) and args is a tuple of the - fixed parameters needed to completely specify the function. + fixed parameters needed to completely specify the function. This + function must be JIT-compiled in `nopython` mode using Numba. bounds: ndarray(float, ndim=2) Sequence of (min, max) pairs for each element in x. diff --git a/quantecon/optimize/tests/__init__.py b/quantecon/optimize/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/quantecon/optimize/tests/test_multivar_max.py b/quantecon/optimize/tests/test_nelder_mead.py similarity index 81% rename from quantecon/optimize/tests/test_multivar_max.py rename to quantecon/optimize/tests/test_nelder_mead.py index c10ac5991..c961a4137 100644 --- a/quantecon/optimize/tests/test_multivar_max.py +++ b/quantecon/optimize/tests/test_nelder_mead.py @@ -8,7 +8,8 @@ from numpy.testing import assert_allclose from nose.tools import raises -from quantecon.optimize import maximize, nelder_mead_algorithm +from quantecon.optimize import nelder_mead +from ..nelder_mead import _nelder_mead_algorithm @njit @@ -38,7 +39,7 @@ def mccormick(x): def bohachevsky(x): # https://www.sfu.ca/~ssurjano/boha.html f = x[0] ** 2 + x[1] ** 2 - 0.3 * np.cos(3 * np.pi * x[0]) - \ - 0.4 * np.cos(4 * np.pi * x[1]) +0.7 + 0.4 * np.cos(4 * np.pi * x[1]) + 0.7 return -f @@ -98,7 +99,7 @@ def colville(x): # https://www.sfu.ca/~ssurjano/colville.html f = 100 * (x[0] ** 2 - x[1]) ** 2 + (x[0] - 1) ** 2 + (x[2] - 1) ** 2 + \ 90 * (x[2] ** 2 - x[3]) ** 2 + 10.1 * \ - ((x[1] - 1)** 2 + (x[3] - 1) ** 2) + 19.8 * (x[1] - 1) * (x[3] - 1) + ((x[1] - 1) ** 2 + (x[3] - 1) ** 2) + 19.8 * (x[1] - 1) * (x[3] - 1) return -f @@ -123,6 +124,7 @@ def goldstein_price(x): f = (1 + p1 * p2) * (30 + p3 * p4) return -f + @njit def sum_squared(x): return - (x ** 2).sum() @@ -153,7 +155,7 @@ def test_rosenbrock(self): x0 = np.array([-2, 1]) - results = maximize(rosenbrock, x0, tol_x=1e-20, tol_f=1e-20) + results = nelder_mead(rosenbrock, x0, tol_x=1e-20, tol_f=1e-20) assert_allclose(results.x, sol, atol=1e-4) assert_allclose(results.fun, fun, atol=1e-4) @@ -164,7 +166,7 @@ def test_powell(self): x0 = np.array([3., -1., 0., 1.]) - results = maximize(powell, x0, tol_x=1e-20, tol_f=1e-20) + results = nelder_mead(powell, x0, tol_x=1e-20, tol_f=1e-20) assert_allclose(results.x, sol, atol=1e-4) assert_allclose(results.fun, fun, atol=1e-4) @@ -177,7 +179,7 @@ def test_mccormick(self): bounds = np.array([[-1.5, 4.], [-3., 4.]]) - results = maximize(mccormick, x0, bounds=bounds) + results = nelder_mead(mccormick, x0, bounds=bounds) assert_allclose(results.x, sol, rtol=1e-3) assert_allclose(results.fun, fun, rtol=1e-3) @@ -189,7 +191,7 @@ def test_bohachevsky(self): # Starting point makes significant difference x0 = np.array([np.pi, -np.pi]) - results = maximize(bohachevsky, x0) + results = nelder_mead(bohachevsky, x0) assert_allclose(results.x, sol, atol=1e-4) assert_allclose(results.fun, fun, atol=1e-4) @@ -200,7 +202,7 @@ def test_easom(self): x0 = np.array([5, -1]) - results = maximize(easom, x0, tol_x=1e-20, tol_f=1e-20) + results = nelder_mead(easom, x0, tol_x=1e-20, tol_f=1e-20) assert_allclose(results.x, sol, atol=1e-4) assert_allclose(results.fun, fun, atol=1e-4) @@ -213,8 +215,8 @@ def test_perm_function(self): sol = np.array([1 / d for d in range(1, int(d)+1)]) fun = 0. - results = maximize(perm_function, x0, bounds=bounds, args=(1., ), - tol_x=1e-30, tol_f=1e-30) + results = nelder_mead(perm_function, x0, bounds=bounds, args=(1., ), + tol_x=1e-30, tol_f=1e-30) assert_allclose(results.x, sol, atol=1e-7) assert_allclose(results.fun, fun, atol=1e-7) @@ -227,8 +229,8 @@ def test_rotated_hyper_ellipsoid(self): sol = np.zeros(d) fun = 0. - results = maximize(rotated_hyper_ellipsoid, x0, bounds=bounds, - tol_x=1e-30, tol_f=1e-30) + results = nelder_mead(rotated_hyper_ellipsoid, x0, bounds=bounds, + tol_x=1e-30, tol_f=1e-30) assert_allclose(results.x, sol, atol=1e-7) assert_allclose(results.fun, fun, atol=1e-7) @@ -239,7 +241,7 @@ def test_booth(self): sol = np.array([1., 3.]) fun = 0. - results = maximize(booth, x0, tol_x=1e-20, tol_f=1e-20) + results = nelder_mead(booth, x0, tol_x=1e-20, tol_f=1e-20) assert_allclose(results.x, sol, atol=1e-7) assert_allclose(results.fun, fun, atol=1e-7) @@ -251,8 +253,8 @@ def test_zakharov(self): sol = np.zeros(5) fun = 0. - results = maximize(zakharov, x0, bounds=bounds, tol_f=1e-30, - tol_x=1e-30) + results = nelder_mead(zakharov, x0, bounds=bounds, tol_f=1e-30, + tol_x=1e-30) assert_allclose(results.x, sol, atol=1e-7) assert_allclose(results.fun, fun, atol=1e-7) @@ -264,8 +266,8 @@ def test_colville(self): sol = np.ones(4) fun = 0. - results = maximize(colville, x0, bounds=bounds, tol_f=1e-35, - tol_x=1e-35) + results = nelder_mead(colville, x0, bounds=bounds, tol_f=1e-35, + tol_x=1e-35) assert_allclose(results.x, sol) assert_allclose(results.fun, fun, atol=1e-7) @@ -278,8 +280,8 @@ def test_styblinski_tang(self): sol = np.array([-2.903534] * d) fun = 39.16599 * d - results = maximize(styblinski_tang, x0, bounds=bounds, tol_f=1e-35, - tol_x=1e-35) + results = nelder_mead(styblinski_tang, x0, bounds=bounds, tol_f=1e-35, + tol_x=1e-35) assert_allclose(results.x, sol, rtol=1e-4) assert_allclose(results.fun, fun, rtol=1e-5) @@ -289,7 +291,7 @@ def test_goldstein_price(self): bounds = np.array([[-2., 2.], [-2., 2.]]) - results = maximize(goldstein_price, x0) + results = nelder_mead(goldstein_price, x0) sol = np.array([0., -1.]) fun = -3. @@ -303,11 +305,10 @@ def test_sum_squared(self): sol = np.zeros(3) fun = 0. - results = maximize(sum_squared, x0, tol_f=1e-50, tol_x=1e-50) + results = nelder_mead(sum_squared, x0, tol_f=1e-50, tol_x=1e-50) assert_allclose(results.x, sol, atol=1e-5) assert_allclose(results.fun, fun, atol=1e-5) - def test_corner_sol(self): sol = np.array([0.]) fun = 0. @@ -315,7 +316,7 @@ def test_corner_sol(self): x0 = np.array([10.]) bounds = np.array([[0., np.inf]]) - results = maximize(f, x0, bounds=bounds, tol_f=1e-20) + results = nelder_mead(f, x0, bounds=bounds, tol_f=1e-20) assert_allclose(results.x, sol) assert_allclose(results.fun, fun) @@ -326,7 +327,7 @@ def test_discontinuous(self): x0 = np.array([-10.]) - results = maximize(g, x0) + results = nelder_mead(g, x0) assert_allclose(results.x, sol) assert_allclose(results.fun, fun) @@ -335,11 +336,11 @@ def test_discontinuous(self): @raises(ValueError) def test_invalid_bounds_1(): vertices = np.array([[-2., 1.], - [1.05 * -2., 1.], - [-2., 1.05 * 1.]]) + [1.05 * -2., 1.], + [-2., 1.05 * 1.]]) x0 = np.array([-2., 1.]) bounds = np.array([[10., -10.], [10., -10.]]) - maximize(rosenbrock, x0, bounds=bounds) + nelder_mead(rosenbrock, x0, bounds=bounds) @raises(ValueError) @@ -349,7 +350,7 @@ def test_invalid_bounds_2(): [-2., 1.05 * 1.]]) x0 = np.array([-2., 1.]) bounds = np.array([[10., -10., 10., -10.]]) - maximize(rosenbrock, x0, bounds=bounds) + nelder_mead(rosenbrock, x0, bounds=bounds) @raises(ValueError) @@ -358,7 +359,7 @@ def test_invalid_ρ(): [1.05 * -2., 1.], [-2., 1.05 * 1.]]) invalid_ρ = -1. - nelder_mead_algorithm(rosenbrock, vertices, ρ=invalid_ρ) + _nelder_mead_algorithm(rosenbrock, vertices, ρ=invalid_ρ) @raises(ValueError) @@ -367,7 +368,7 @@ def test_invalid_χ(): [1.05 * -2., 1.], [-2., 1.05 * 1.]]) invalid_χ = 0.5 - nelder_mead_algorithm(rosenbrock, vertices, χ=invalid_χ) + _nelder_mead_algorithm(rosenbrock, vertices, χ=invalid_χ) @raises(ValueError) @@ -377,7 +378,7 @@ def test_invalid_ρχ(): [-2., 1.05 * 1.]]) ρ = 2 χ = 1.5 - nelder_mead_algorithm(rosenbrock, vertices, ρ=ρ, χ=χ) + _nelder_mead_algorithm(rosenbrock, vertices, ρ=ρ, χ=χ) @raises(ValueError) @@ -386,7 +387,7 @@ def test_invalid_γ(): [1.05 * -2., 1.], [-2., 1.05 * 1.]]) invalid_γ = -1e-7 - nelder_mead_algorithm(rosenbrock, vertices, γ=invalid_γ) + _nelder_mead_algorithm(rosenbrock, vertices, γ=invalid_γ) @raises(ValueError) @@ -395,7 +396,7 @@ def test_invalid_σ(): [1.05 * -2., 1.], [-2., 1.05 * 1.]]) invalid_σ = 1. + 1e-7 - nelder_mead_algorithm(rosenbrock, vertices, σ=invalid_σ) + _nelder_mead_algorithm(rosenbrock, vertices, σ=invalid_σ) if __name__ == '__main__':