From 852e21b1e72a1a27b69bef72524190601f345564 Mon Sep 17 00:00:00 2001 From: Damar Wicaksono Date: Tue, 10 Dec 2024 13:38:27 +0100 Subject: [PATCH] Add the Rosenbrock function for optimization and metamodeling The M-dimensional Rosenbrock function has been added to the codebase. The function is a widely used function for testing optimization algorithms. There is also a known use case in a metamodeling exercise. The documentation has been updated accordingly. This commit should resolve Issue #405. --- CHANGELOG.md | 2 + docs/_toc.yml | 2 + docs/fundamentals/metamodeling.md | 1 + docs/fundamentals/optimization.md | 9 +- docs/references.bib | 33 ++ docs/test-functions/available.md | 1 + docs/test-functions/rosenbrock.md | 302 ++++++++++++++++++ src/uqtestfuns/test_functions/__init__.py | 2 + src/uqtestfuns/test_functions/rosenbrock.py | 155 +++++++++ .../builtin_test_functions/test_rosenbrock.py | 55 ++++ 10 files changed, 558 insertions(+), 4 deletions(-) create mode 100644 docs/test-functions/rosenbrock.md create mode 100644 src/uqtestfuns/test_functions/rosenbrock.py create mode 100644 tests/builtin_test_functions/test_rosenbrock.py diff --git a/CHANGELOG.md b/CHANGELOG.md index a7f70a9..ea8ae6f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- The M-dimensional Rosenbrock function for optimization and metamodeling + exercises. - An additional use case reference for the Ackley function (as a metamodeling test function). - The ten-dimensional inert function from Linkletter et al. (2006) with diff --git a/docs/_toc.yml b/docs/_toc.yml index 758e83d..01cc95a 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -158,6 +158,8 @@ parts: title: Simple Portfolio Model - file: test-functions/robot-arm title: Robot Arm + - file: test-functions/rosenbrock + title: Rosenbrock - file: test-functions/rs-circular-bar title: RS - Circular Bar - file: test-functions/rs-quadratic diff --git a/docs/fundamentals/metamodeling.md b/docs/fundamentals/metamodeling.md index 09ce403..588270f 100644 --- a/docs/fundamentals/metamodeling.md +++ b/docs/fundamentals/metamodeling.md @@ -60,6 +60,7 @@ in the comparison of metamodeling approaches. | {ref}`OTL Circuit ` | 6 / 20 | `OTLCircuit()` | | {ref}`Piston Simulation ` | 7 / 20 | `Piston()` | | {ref}`Robot Arm ` | 8 | `RobotArm()` | +| {ref}`Rosenbrock ` | M | `Rosenbrock()` | | {ref}`Solar Cell Model ` | 5 | `SolarCell()` | | {ref}`Sulfur ` | 9 | `Sulfur()` | | {ref}`Undamped Oscillator ` | 6 | `UndampedOscillator()` | diff --git a/docs/fundamentals/optimization.md b/docs/fundamentals/optimization.md index fb3a18a..16ca72f 100644 --- a/docs/fundamentals/optimization.md +++ b/docs/fundamentals/optimization.md @@ -18,10 +18,11 @@ kernelspec: The table below listed the available test functions typically used in the comparison of global optimization methods. -| Name | Input Dimension | Constructor | -|:-----------------------------------------------------------:|:---------------:|:--------------------:| -| {ref}`Ackley ` | M | `Ackley()` | -| {ref}`Forrester et al. (2008) ` | 1 | `Forrester2008()` | +| Name | Input Dimension | Constructor | +|:---------------------------------------------------------:|:---------------:|:-----------------:| +| {ref}`Ackley ` | M | `Ackley()` | +| {ref}`Forrester et al. (2008) ` | 1 | `Forrester2008()` | +| {ref}`Rosenbrock ` | M | `Rosenbrock()` | In a Python terminal, you can list all the available functions relevant for optimization applications using ``list_functions()`` and filter the results diff --git a/docs/references.bib b/docs/references.bib index d279184..6f7f404 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -1065,4 +1065,37 @@ @Article{Kaintura2017 doi = {10.1007/s00366-017-0507-0}, } +@Article{Rosenbrock1960, + author = {Rosenbrock, H. H.}, + journal = {The Computer Journal}, + title = {An automatic method for finding the greatest or least value of a function}, + year = {1960}, + number = {3}, + pages = {175--184}, + volume = {3}, + doi = {10.1093/comjnl/3.3.175}, +} + +@Article{Picheny2013, + author = {Picheny, Victor and Wagner, Tobias and Ginsbourger, David}, + journal = {Structural and Multidisciplinary Optimization}, + title = {A benchmark of kriging-based infill criteria for noisy optimization}, + year = {2013}, + number = {3}, + pages = {607--626}, + volume = {48}, + doi = {10.1007/s00158-013-0919-4}, +} + +@Article{Tan2015, + author = {Tan, Matthias Hwai Yong}, + journal = {Technometrics}, + title = {Stochastic polynomial interpolation for uncertainty quantification with computer experiments}, + year = {2015}, + number = {4}, + pages = {457--467}, + volume = {57}, + doi = {10.1080/00401706.2014.950431}, +} + @Comment{jabref-meta: databaseType:bibtex;} diff --git a/docs/test-functions/available.md b/docs/test-functions/available.md index adb65b6..182ec95 100644 --- a/docs/test-functions/available.md +++ b/docs/test-functions/available.md @@ -82,6 +82,7 @@ regardless of their typical applications. | {ref}`Piston Simulation ` | 7 / 20 | `Piston()` | | {ref}`Simple Portfolio Model ` | 3 | `Portfolio3D()` | | {ref}`Robot Arm ` | 8 | `RobotArm()` | +| {ref}`Rosenbrock ` | M | `Rosenbrock()` | | {ref}`RS - Circular Bar ` | 2 | `RSCircularBar()` | | {ref}`RS - Quadratic ` | 2 | `RSQuadratic()` | | {ref}`SaltelliLinear ` | M | `SaltelliLinear()` | diff --git a/docs/test-functions/rosenbrock.md b/docs/test-functions/rosenbrock.md new file mode 100644 index 0000000..97282b5 --- /dev/null +++ b/docs/test-functions/rosenbrock.md @@ -0,0 +1,302 @@ +--- +jupytext: + formats: ipynb,md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.1 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +(test-functions:rosenbrock)= +# Rosenbrock Function + +The Rosenbrock function, originally introduced in {cite}`Rosenbrock1960` +as a two-dimensional scalar-valued test function for global optimization, +was later generalized to $M$ dimensions. +It has since become a widely used benchmark for global optimization methods +(e.g., {cite}`Dixon1978, Picheny2013`). +In {cite}`Tan2015`, the function was employed in a metamodeling exercise. + +The function is also known as the valley function or the banana function. + +```{code-cell} ipython3 +import numpy as np +import matplotlib.pyplot as plt +import uqtestfuns as uqtf +``` + +The surface and contour plots for the two-dimensional Rosenbrock function are +shown below for the default parameter set and for $x \in [-2, 2] \times [-1, 3]$. + +As shown, the function features a curved, non-convex valley. +While it is relatively easy to reach the valley, the convergence to the global +minimum is difficult. + +```{code-cell} ipython3 +:tags: [remove-input] + +from mpl_toolkits.axes_grid1 import make_axes_locatable + +# --- Create 2D data +x1_1d = np.arange(-2, 2, 0.05) +x2_1d = np.arange(-1, 3, 0.05) +my_fun = uqtf.Rosenbrock(input_dimension=2) +mesh_2d = np.meshgrid(x1_1d, x2_1d) +xx_2d = np.array(mesh_2d).T.reshape(-1, 2) +yy_2d = my_fun(xx_2d) + +# --- Create two-dimensional plots +fig = plt.figure(figsize=(10, 10)) + +# Surface +axs_1 = plt.subplot(121, projection='3d') +axs_1.plot_surface( + mesh_2d[0], + mesh_2d[1], + yy_2d.reshape(len(x1_1d), len(x2_1d)).T, + linewidth=0, + cmap="plasma", + antialiased=False, + alpha=0.5 +) +axs_1.set_xlabel("$x_1$", fontsize=14) +axs_1.set_ylabel("$x_2$", fontsize=14) +axs_1.set_title("Surface plot of 2D Rosenbrock", fontsize=14) + +# Contour +axs_2 = plt.subplot(122) +cf = axs_2.contour( + mesh_2d[0], mesh_2d[1], yy_2d.reshape(len(x1_1d), len(x2_1d)).T, + levels=1000, cmap="plasma", +) +axs_2.set_xlabel("$x_1$", fontsize=14) +axs_2.set_ylabel("$x_2$", fontsize=14) +axs_2.set_title("Contour plot of 2D Rosenbrock", fontsize=14) +divider = make_axes_locatable(axs_2) +cax = divider.append_axes('right', size='5%', pad=0.05) +fig.colorbar(cf, cax=cax, orientation='vertical') +axs_2.axis('scaled') + +fig.tight_layout(pad=3.0) +plt.gcf().set_dpi(75); +``` + +## Test function instance + +To create a default instance of the test function, type: + +```{code-cell} ipython3 +my_testfun = uqtf.Rosenbrock() +``` + +Check if it has been correctly instantiated: + +```{code-cell} ipython3 +print(my_testfun) +``` + +By default, the input dimension is set to $2$[^default_dimension]. +To create an instance with another value of input dimension, +pass an integer to the parameter `input_dimension` (keyword only). +For example, to create an instance of 10-dimensional Rosenbrock function, type: + +```python +my_testfun = uqtf.Ackley(input_dimension=10) +``` + +## Description + +The generalized Rosenbrock function is defined as follows: + +$$ +\mathcal{M}(\boldsymbol{x}; \boldsymbol{p}) = \frac{1}{d} \left[ +\left( \sum_{i = 1}^{M - 1} (x_i - a)^2 + b \, (x_{i + 1} - x_i^2)^2 \right) +- c \right] +$$ + +where $\boldsymbol{x} = \{ x_1, \ldots, x_M \}$ is the $M$-dimensional vector +of input variables, as defined below, and +$\boldsymbol{p} = \{ a, b, c, d \}$ is the set of parameters also defined +further below. + +```{info} +For $M < 2$, the Rosenbrock function returns a constant $0$ for any $boldsymbol{x}$. +``` + +## Input + +The Rosenbrock function is defined on $\mathbb{R}^M$, but in the context +of optimization test function, the search space is often limited as shown +in the table below. + +```{code-cell} ipython3 +:tags: [hide-input] + +print(my_testfun.prob_input) +``` + +## Parameters + +The Rosenbrock function requires four additional parameters +to complete the specification. +The default values are shown below. + +```{code-cell} ipython3 +:tags: [hide-input] + +print(my_testfun.parameters) +``` + +The parameter $a$ controls the location and value of the global optimum. +The parameter $b$ controls the steepness and width of the valley. +In particular, it determines the scale of variation of the function; large +value of $b$ creates a steeper and tight valley. + +The parameters $c$ and $d$ are used to shift and scale the function such that +its mean and standard deviation become $0.0$ and $1.0$, respectively. + +Below are some contour plots of the function in two dimensions with different +values of $a$ and $b$ along with the location of the global optimum. + +```{code-cell} ipython3 +:tags: [remove-input] + +from mpl_toolkits.axes_grid1 import make_axes_locatable + +# --- Create 2D data +x1_1d = np.arange(-2, 2, 0.05) +x2_1d = np.arange(-1, 3, 0.05) +mesh_2d = np.meshgrid(x1_1d, x2_1d) +xx_2d = np.array(mesh_2d).T.reshape(-1, 2) + +# --- Create contour plots +fig, axs = plt.subplots(2, 2, figsize=(10, 10)) + +# a = 1.0, b = 100.0 +my_fun = uqtf.Rosenbrock(input_dimension=2) +yy_2d = my_fun(xx_2d) +a = my_fun.parameters["a"] +b = my_fun.parameters["b"] + +cf = axs[0, 0].contour( + mesh_2d[0], mesh_2d[1], yy_2d.reshape(len(x1_1d), len(x2_1d)).T, + levels=1000, cmap="plasma", +) +axs[0, 0].set_xlabel("$x_1$", fontsize=14) +axs[0, 0].set_ylabel("$x_2$", fontsize=14) +axs[0, 0].set_title(f"a = {a}, b = {b}", fontsize=14) +divider = make_axes_locatable(axs[0, 0]) +cax = divider.append_axes('right', size='5%', pad=0.05) +fig.colorbar(cf, cax=cax, orientation='vertical') +axs[0, 0].axis('scaled') + +axs[0, 0].scatter(a, a**2, marker="x", s=150, color="k") + +# a = 1.6, b = 100.0 +my_fun = uqtf.Rosenbrock(input_dimension=2) +my_fun.parameters["a"] = 1.6 +yy_2d = my_fun(xx_2d) +a = my_fun.parameters["a"] +b = my_fun.parameters["b"] + +cf = axs[0, 1].contour( + mesh_2d[0], mesh_2d[1], yy_2d.reshape(len(x1_1d), len(x2_1d)).T, + levels=1000, cmap="plasma", +) +axs[0, 1].set_xlabel("$x_1$", fontsize=14) +axs[0, 1].set_ylabel("$x_2$", fontsize=14) +axs[0, 1].set_title(f"a = {a}, b = {b}", fontsize=14) +divider = make_axes_locatable(axs[0, 1]) +cax = divider.append_axes('right', size='5%', pad=0.05) +fig.colorbar(cf, cax=cax, orientation='vertical') +axs[0, 1].axis('scaled') + +axs[0, 1].scatter(a, a**2, marker="x", s=150, color="k") + +# a = 1.0, b = 500.0 +my_fun = uqtf.Rosenbrock(input_dimension=2) +my_fun.parameters["b"] = 500.0 +yy_2d = my_fun(xx_2d) +a = my_fun.parameters["a"] +b = my_fun.parameters["b"] + +cf = axs[1, 0].contour( + mesh_2d[0], mesh_2d[1], yy_2d.reshape(len(x1_1d), len(x2_1d)).T, + levels=1000, cmap="plasma", +) +axs[1, 0].set_xlabel("$x_1$", fontsize=14) +axs[1, 0].set_ylabel("$x_2$", fontsize=14) +axs[1, 0].set_title(f"a = {a}, b = {b}", fontsize=14) +divider = make_axes_locatable(axs[1, 0]) +cax = divider.append_axes('right', size='5%', pad=0.05) +fig.colorbar(cf, cax=cax, orientation='vertical') +axs[1, 0].axis('scaled') + +axs[1, 0].scatter(a, a**2, marker="x", s=150, color="k") + +# a = 1.6, b = 500.0 +my_fun = uqtf.Rosenbrock(input_dimension=2) +my_fun.parameters["a"] = 1.6 +my_fun.parameters["b"] = 500.0 +yy_2d = my_fun(xx_2d) +a = my_fun.parameters["a"] +b = my_fun.parameters["b"] + +cf = axs[1, 1].contour( + mesh_2d[0], mesh_2d[1], yy_2d.reshape(len(x1_1d), len(x2_1d)).T, + levels=1000, cmap="plasma", +) +axs[1, 1].set_xlabel("$x_1$", fontsize=14) +axs[1, 1].set_ylabel("$x_2$", fontsize=14) +axs[1, 1].set_title(f"a = {a}, b = {b}", fontsize=14) +divider = make_axes_locatable(axs[1, 1]) +cax = divider.append_axes('right', size='5%', pad=0.05) +fig.colorbar(cf, cax=cax, orientation='vertical') +axs[1, 1].axis('scaled') + +axs[1, 1].scatter(a, a**2, marker="x", s=150, color="k") + +fig.tight_layout() +plt.gcf().set_dpi(75); +``` + +As mentioned earlier, the changing the value of $a$ modifies the location +of the global optimum (and in $M > 2$ modifies the function value at the +optimum). Changing the value of $b$, on the other hand, modifies the scale +of the function variation and the steepness of the valley (see the color bar). + +## Reference results + +This section provides several reference results related to the test function. + +### Optimum values + +The global optimum of the Rosenbrock function is located at + +$$ +\begin{aligned} +\boldsymbol{x}^* & = \left( x^*_1, \ldots, x^*_M \right), \\ +x_i^* & = a^{2^{i - 1}}. +\end{aligned} +$$ + +In dimension two, the function value at the optimum location is always $0.0$ +(but not in higher dimension!). + +## References + +```{bibliography} +:style: unsrtalpha +:filter: docname in docnames +``` + +[^default_dimension]: This default dimension applies to all variable dimension +test functions. It will be used if the `input_dimension` argument is not given. + +[^location]: The original two-dimensional expression can be found in Eq. (10), +Section 3.2, in {cite}`Rosenbrock1960`. diff --git a/src/uqtestfuns/test_functions/__init__.py b/src/uqtestfuns/test_functions/__init__.py index 39e0c71..4f5f06e 100644 --- a/src/uqtestfuns/test_functions/__init__.py +++ b/src/uqtestfuns/test_functions/__init__.py @@ -49,6 +49,7 @@ from .piston import Piston from .portfolio_3d import Portfolio3D from .robot_arm import RobotArm +from .rosenbrock import Rosenbrock from .rs_circular_bar import RSCircularBar from .rs_quadratic import RSQuadratic from .saltelli_linear import SaltelliLinear @@ -126,6 +127,7 @@ "Piston", "Portfolio3D", "RobotArm", + "Rosenbrock", "RSCircularBar", "RSQuadratic", "SaltelliLinear", diff --git a/src/uqtestfuns/test_functions/rosenbrock.py b/src/uqtestfuns/test_functions/rosenbrock.py new file mode 100644 index 0000000..aab46be --- /dev/null +++ b/src/uqtestfuns/test_functions/rosenbrock.py @@ -0,0 +1,155 @@ +""" +This module provides an implementation of the Rosenbrock function. + +The Rosenbrock function, originally introduced in [1] as a two-dimensional +scalar-valued test function for global optimization, was later generalized +to M dimensions. +It has since become a widely used benchmark for global optimization methods +(e.g., [2], [3]). +In [4], the function was employed in a metamodeling exercise. + +The function features a curved, non-convex valley. +While it is relatively easy to reach the valley, the convergence to the global +minimum is difficult. + +References +---------- +1. H. H. Rosenbrock, “An Automatic Method for Finding the Greatest or Least + Value of a Function,” The Computer Journal, vol. 3, no. 3, pp. 175–184, + 1960. + DOI: 10.1093/comjnl/3.3.175 +2. L. C. W. Dixon and G. P. Szegö. Towards global optimization 2, + chapter The global optimization problem: an introduction, pages 1–15. + North-Holland, Amsterdam, 1978. +3. V. Picheny, T. Wagner, and D. Ginsbourger, “A benchmark of kriging-based + infill criteria for noisy optimization,” Structural and Multidisciplinary + Optimization, vol. 48, no. 3, pp. 607–626, 2013. + DOI: 10.1007/s00158-013-0919-4 +4. M. H. Y. Tan, “Stochastic Polynomial Interpolation for Uncertainty + Quantification With Computer Experiments,” Technometrics, vol. 57, no. 4, + pp. 457–467, 2015. + DOI: 10.1080/00401706.2014.950431 +""" + +import numpy as np + +from uqtestfuns.core.uqtestfun_abc import UQTestFunVarDimABC +from uqtestfuns.core.custom_typing import ProbInputSpecs, FunParamSpecs + +__all__ = ["Rosenbrock"] + + +AVAILABLE_INPUTS: ProbInputSpecs = { + "Picheny2013": { + "function_id": "Rosenbrock", + "description": ( + "Search domain for the Rosenbrock function from Picheny et al. " + "(2013)" + ), + "marginals": [ + { + "name": "x", + "distribution": "uniform", + "parameters": [-5.0, 10.0], + "description": None, + }, + ], + "copulas": None, + } +} + + +AVAILABLE_PARAMETERS: FunParamSpecs = { + "Rosenbrock1960": { + "function_id": "Rosenbrock", + "description": ( + "Parameter set for the Rosenbrock function from Rosenbrock (1960)" + ), + "declared_parameters": [ + { + "keyword": "a", + "value": 1.0, + "type": float, + "description": "Global optimum location and value", + }, + { + "keyword": "b", + "value": 100.0, + "type": float, + "description": ( + "Steepness, curvature, and width of the valley" + ), + }, + { + "keyword": "c", + "value": 0.0, + "type": float, + "description": "Shift parameter", + }, + { + "keyword": "d", + "value": 1.0, + "type": float, + "description": "Scale parameter", + }, + ], + }, +} + + +def evaluate( + xx: np.ndarray, a: float, b: float, c: float, d: float +) -> np.ndarray: + """Evaluate the Rosenbrock function on a set of input values. + + Parameters + ---------- + xx : np.ndarray + M-Dimensional input values given by an N-by-M array where + N is the number of input values. + a : float + The parameter controlling the magnitude of the quadratic penalty term + between a given input and the parameter; the global optimum location + and value are controlled by this parameter. + b : float + The parameter controlling the magnitude of the quadratic penalty term + between two adjacent input variables; the steepness, curvature, and + width of the valley are controlled by this parameter. + c : float + The shift parameter. + d : float + The scale parameter. + + Returns + ------- + np.ndarray + The output of the test function evaluated on the input values. + The output is a 1-dimensional array of length N. + + Notes + ----- + - The function returns a constant zero for M < 2. + """ + if xx.shape[1] == 1: + return np.zeros(xx.shape[0]) + + yy = (xx[:, :-1] - a) ** 2 + b * (xx[:, 1:] - xx[:, :-1] ** 2) ** 2 + yy = np.sum(yy, axis=1) + # Shift and rescale + yy = (yy - c) / d + + return yy + + +class Rosenbrock(UQTestFunVarDimABC): + """A concrete implementation of the Rosenbrock test function.""" + + _tags = ["optimization", "metamodeling"] + _description = ( + "Optimization test function from Rosenbrock (1960), " + "also known as the banana function" + ) + _available_inputs = AVAILABLE_INPUTS + _available_parameters = AVAILABLE_PARAMETERS + + evaluate = staticmethod(evaluate) # type: ignore diff --git a/tests/builtin_test_functions/test_rosenbrock.py b/tests/builtin_test_functions/test_rosenbrock.py new file mode 100644 index 0000000..a8247dd --- /dev/null +++ b/tests/builtin_test_functions/test_rosenbrock.py @@ -0,0 +1,55 @@ +""" +Test module for the Ackley test function. + +Notes +----- +- The tests defined in this module deals with + the correctness of the evaluation. +""" + +import numpy as np +import pytest + +from uqtestfuns import Rosenbrock + + +def test_one_dimensional(): + """Test evaluating the one-dimensional Rosenbrock function.""" + my_fun = Rosenbrock(input_dimension=1) + + xx = my_fun.prob_input.get_sample(100) + yy = my_fun(xx) + + # One-dimensional Rosenbrock always returns zero + assert np.allclose(yy, 0.0) + + +@pytest.mark.parametrize("input_dimension", [2, 3, 10]) +def test_optimum_value_a_eq_1(input_dimension): + """Test the optimum value, regardless of the dimension, when a == 1.""" + my_fun = Rosenbrock(input_dimension=input_dimension) + + xx = np.ones((1, my_fun.input_dimension)) + yy = my_fun(xx) + + # For a == 0, The optima of the Ackley function is at 1.0s + assert np.allclose(yy, 0.0) + + +@pytest.mark.parametrize("input_dimension", [2, 3, 10]) +def test_optimum_value_a_eq_0(input_dimension): + """Test the optimum value, regardless of the dimension, when a == 0""" + my_fun = Rosenbrock(input_dimension=input_dimension) + my_fun.parameters["a"] = 0.0 + + # The optima of the function when a = 0 is not at 1.0's + xx = np.ones((1, my_fun.input_dimension)) + yy = my_fun(xx) + + assert not np.allclose(yy, 0.0) + + # but at 0.0's + xx = np.zeros((1, my_fun.input_dimension)) + yy = my_fun(xx) + + assert np.allclose(yy, 0.0)