Nelder-Mead Generator adapted from SciPy¶

Most of the algorithms in scipy.optimize are self-contained functions that operate on the user-provided func. Xopt has adapted the Nelder-Mead directly from scipy.optimize to be in a generator form. This allows for the manual stepping through the algorithm.

In [1]:

from xopt.generators.scipy.neldermead import NelderMeadGenerator
from xopt import Evaluator, VOCS
from xopt.resources.test_functions.rosenbrock import rosenbrock

import pandas as pd

from xopt import Xopt
import numpy as np

from scipy.optimize import fmin

# from xopt import output_notebook
# output_notebook()

import matplotlib.pyplot as plt

from xopt.generators.scipy.neldermead import NelderMeadGenerator from xopt import Evaluator, VOCS from xopt.resources.test_functions.rosenbrock import rosenbrock import pandas as pd from xopt import Xopt import numpy as np from scipy.optimize import fmin # from xopt import output_notebook # output_notebook() import matplotlib.pyplot as plt

Nelder-Mead optimization of the Rosenbrock function with Xopt¶

In [2]:

YAML = """
generator:
  name: neldermead
  initial_point: {x0: -1, x1: -1}
  adaptive: true
  xatol: 0.0001
  fatol: 0.0001  
evaluator:
  function: xopt.resources.test_functions.rosenbrock.evaluate_rosenbrock
vocs:
  variables:
    x0: [-5, 5]
    x1: [-5, 5]
  objectives: {y: MINIMIZE}
"""
X = Xopt.from_yaml(YAML)

YAML = """ generator: name: neldermead initial_point: {x0: -1, x1: -1} adaptive: true xatol: 0.0001 fatol: 0.0001 evaluator: function: xopt.resources.test_functions.rosenbrock.evaluate_rosenbrock vocs: variables: x0: [-5, 5] x1: [-5, 5] objectives: {y: MINIMIZE} """ X = Xopt.from_yaml(YAML)

In [3]:

XMIN = [1, 1]  # True minimum

XMIN = [1, 1] # True minimum

In [4]:

X.run()
X.data

X.run() X.data

Out[4]:

	x0	x1	y	xopt_runtime	xopt_error
0	-1.000000	-1.000000	4.040000e+02	0.000008	False
1	-1.050000	-1.000000	4.462531e+02	0.000007	False
2	-1.000000	-1.050000	4.242500e+02	0.000005	False
3	-0.950000	-1.050000	3.850281e+02	0.000005	False
4	-0.900000	-1.075000	3.589325e+02	0.000006	False
...	...	...	...	...	...
120	0.999935	0.999867	5.114951e-09	0.000005	False
121	0.999877	0.999764	2.587916e-08	0.000005	False
122	0.999999	0.999995	5.309344e-10	0.000005	False
123	1.000045	1.000097	7.751675e-09	0.000005	False
124	0.999963	0.999925	1.412126e-09	0.000005	False

125 rows × 5 columns

In [5]:

# Evaluation progression
X.data["y"].plot(marker=".")
plt.yscale("log")
plt.xlabel("iteration")
plt.ylabel("Rosenbrock value")

# Evaluation progression X.data["y"].plot(marker=".") plt.yscale("log") plt.xlabel("iteration") plt.ylabel("Rosenbrock value")

Out[5]:

Text(0, 0.5, 'Rosenbrock value')

In [6]:

# Minimum
dict(X.data.iloc[X.data["y"].argmin()])

# Minimum dict(X.data.iloc[X.data["y"].argmin()])

Out[6]:

{'x0': np.float64(0.9999988592114838),
 'x1': np.float64(0.9999954170486077),
 'y': np.float64(5.309343918637161e-10),
 'xopt_runtime': np.float64(5.259999852569308e-06),
 'xopt_error': np.False_}

Visualize¶

In [7]:

fig, ax = plt.subplots(figsize=(8, 8))

Xgrid, Ygrid = np.meshgrid(np.linspace(-2, 2, 201), np.linspace(-2, 2, 201))

Zgrid = np.vectorize(lambda x, y: rosenbrock([x, y]))(Xgrid, Ygrid)
Zgrid = np.log(Zgrid + 1)

ax.pcolormesh(Xgrid, Ygrid, Zgrid)
ax.contour(Xgrid, Ygrid, Zgrid, levels=10, colors="black")
ax.set_xlabel("x0")
ax.set_ylabel("x1")


# Add all evaluations
ax.plot(X.data["x0"], X.data["x1"], color="red", alpha=0.5, marker=".")
ax.scatter(XMIN[0], XMIN[1], 50, marker="o", color="orange", label="True minimum")
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
# plt.legend()
ax.set_title("Xopt's Nelder-Mead progression")

fig, ax = plt.subplots(figsize=(8, 8)) Xgrid, Ygrid = np.meshgrid(np.linspace(-2, 2, 201), np.linspace(-2, 2, 201)) Zgrid = np.vectorize(lambda x, y: rosenbrock([x, y]))(Xgrid, Ygrid) Zgrid = np.log(Zgrid + 1) ax.pcolormesh(Xgrid, Ygrid, Zgrid) ax.contour(Xgrid, Ygrid, Zgrid, levels=10, colors="black") ax.set_xlabel("x0") ax.set_ylabel("x1") # Add all evaluations ax.plot(X.data["x0"], X.data["x1"], color="red", alpha=0.5, marker=".") ax.scatter(XMIN[0], XMIN[1], 50, marker="o", color="orange", label="True minimum") ax.set_xlim(-2, 2) ax.set_ylim(-2, 2) # plt.legend() ax.set_title("Xopt's Nelder-Mead progression")

Out[7]:

Text(0.5, 1.0, "Xopt's Nelder-Mead progression")

In [8]:

# Manually step the algorithm and collect simplexes
X = Xopt.from_yaml(YAML)
simplexes = []
while not X.generator.is_done:
    X.step()
    simplexes.append(X.generator.simplex)

# Manually step the algorithm and collect simplexes X = Xopt.from_yaml(YAML) simplexes = [] while not X.generator.is_done: X.step() simplexes.append(X.generator.simplex)

In [9]:

def plot_simplex(simplex, ax=None):
    x0 = simplex["x0"]
    x1 = simplex["x1"]
    x0 = np.append(x0, x0[0])
    x1 = np.append(x1, x1[0])
    ax.plot(x0, x1)


fig, ax = plt.subplots(figsize=(8, 8))
ax.pcolormesh(Xgrid, Ygrid, Zgrid)
# ax.contour(Xgrid, Ygrid, Zgrid, levels=10, colors='black')
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_xlabel("x0")
ax.set_ylabel("x1")
ax.set_title("Nelder-Mead simplex progression")

ax.scatter(XMIN[0], XMIN[1], 50, marker="o", color="orange", label="True minimum")

for simplex in simplexes:
    plot_simplex(simplex, ax)

def plot_simplex(simplex, ax=None): x0 = simplex["x0"] x1 = simplex["x1"] x0 = np.append(x0, x0[0]) x1 = np.append(x1, x1[0]) ax.plot(x0, x1) fig, ax = plt.subplots(figsize=(8, 8)) ax.pcolormesh(Xgrid, Ygrid, Zgrid) # ax.contour(Xgrid, Ygrid, Zgrid, levels=10, colors='black') ax.set_xlim(-2, 2) ax.set_ylim(-2, 2) ax.set_xlabel("x0") ax.set_ylabel("x1") ax.set_title("Nelder-Mead simplex progression") ax.scatter(XMIN[0], XMIN[1], 50, marker="o", color="orange", label="True minimum") for simplex in simplexes: plot_simplex(simplex, ax)

Compare with scipy.optimize.fmin Nelder-Mead¶

Notice that fmin is much faster here. This is because the function runs very fast, so the internal Xopt bookkeeping overhead dominates.

In [10]:

result = fmin(rosenbrock, [-1, -1])
result

result = fmin(rosenbrock, [-1, -1]) result

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 67
         Function evaluations: 125

Out[10]:

array([0.99999886, 0.99999542])

In [11]:

X = Xopt.from_yaml(YAML)

X = Xopt.from_yaml(YAML)

In [12]:

X.run()
# Almost exactly the same number evaluations.
len(X.data)

X.run() # Almost exactly the same number evaluations. len(X.data)

Out[12]:

125

In [13]:

# results are the same
xbest = X.data.iloc[X.data["y"].argmin()]
xbest["x0"] == result[0], xbest["x1"] == result[1]

# results are the same xbest = X.data.iloc[X.data["y"].argmin()] xbest["x0"] == result[0], xbest["x1"] == result[1]

Out[13]:

(np.True_, np.True_)

NelderMeadGenerator object¶

In [14]:

NelderMeadGenerator.model_fields

NelderMeadGenerator.model_fields

Out[14]:

{'supports_batch_generation': FieldInfo(annotation=bool, required=False, default=False, description='flag that describes if this generator can generate batches of points', exclude=True, frozen=True),
 'supports_multi_objective': FieldInfo(annotation=bool, required=False, default=False, description='flag that describes if this generator can solve multi-objective problems', exclude=True, frozen=True),
 'vocs': FieldInfo(annotation=VOCS, required=True, description='generator VOCS', exclude=True),
 'data': FieldInfo(annotation=Union[DataFrame, NoneType], required=False, default=None, description='generator data', exclude=True),
 'initial_point': FieldInfo(annotation=Union[Dict[str, float], NoneType], required=False, default=None),
 'initial_simplex': FieldInfo(annotation=Union[Dict[str, Union[List[float], ndarray]], NoneType], required=False, default=None),
 'adaptive': FieldInfo(annotation=bool, required=False, default=True, description='Change hyperparameters based on dimensionality'),
 'xatol': FieldInfo(annotation=float, required=False, default=0.0001, description='Tolerance in x value'),
 'fatol': FieldInfo(annotation=float, required=False, default=0.0001, description='Tolerance in function value'),
 'current_state': FieldInfo(annotation=SimplexState, required=False, default=SimplexState(astg=-1, N=None, kend=0, jend=0, ind=None, sim=None, fsim=None, fxr=None, x=None, xr=None, xe=None, xc=None, xcc=None, xbar=None, doshrink=0, ngen=0)),
 'future_state': FieldInfo(annotation=Union[SimplexState, NoneType], required=False, default=None),
 'x': FieldInfo(annotation=Union[ndarray, NoneType], required=False, default=None),
 'y': FieldInfo(annotation=Union[float, NoneType], required=False, default=None),
 'is_done_bool': FieldInfo(annotation=bool, required=False, default=False)}

In [15]:

Xbest = [33, 44]


def f(inputs, verbose=False):
    if verbose:
        print(f"evaluate f({inputs})")
    x0 = inputs["x0"]
    x1 = inputs["x1"]

    # if x0 < 10:
    #    raise ValueError('test XXXX')

    y = (x0 - Xbest[0]) ** 2 + (x1 - Xbest[1]) ** 2

    return {"y": y}


ev = Evaluator(function=f)
vocs = VOCS(
    variables={"x0": [-100, 100], "x1": [-100, 100]}, objectives={"y": "MINIMIZE"}
)
vocs.json()

Xbest = [33, 44] def f(inputs, verbose=False): if verbose: print(f"evaluate f({inputs})") x0 = inputs["x0"] x1 = inputs["x1"] # if x0 < 10: # raise ValueError('test XXXX') y = (x0 - Xbest[0]) ** 2 + (x1 - Xbest[1]) ** 2 return {"y": y} ev = Evaluator(function=f) vocs = VOCS( variables={"x0": [-100, 100], "x1": [-100, 100]}, objectives={"y": "MINIMIZE"} ) vocs.json()

Out[15]:

'{"variables":{"x0":[-100.0,100.0],"x1":[-100.0,100.0]},"constraints":{},"objectives":{"y":"MINIMIZE"},"constants":{},"observables":[]}'

In [16]:

# check output
f(vocs.random_inputs()[0])

# check output f(vocs.random_inputs()[0])

Out[16]:

{'y': 13023.948448483678}

In [17]:

G = NelderMeadGenerator(vocs=vocs)
inputs = G.generate(1)
inputs

G = NelderMeadGenerator(vocs=vocs) inputs = G.generate(1) inputs

Out[17]:

[{'x0': np.float64(-27.020540092972254), 'x1': np.float64(29.869929860151046)}]

In [18]:

# Further generate calls will continue to produce same point, as with BO
G.generate(1)

# Further generate calls will continue to produce same point, as with BO G.generate(1)

Out[18]:

[{'x0': np.float64(-27.020540092972254), 'x1': np.float64(29.869929860151046)}]

In [19]:

ev.evaluate(inputs[0])

ev.evaluate(inputs[0])

Out[19]:

{'y': np.float64(3802.1241152091407),
 'xopt_runtime': 4.057999831275083e-06,
 'xopt_error': False}

In [20]:

# Adding new data will advance state to next step, and next generate() will yield new point
G.add_data(pd.DataFrame([ev.evaluate(inputs[0])]))
G.generate(1)

# Adding new data will advance state to next step, and next generate() will yield new point G.add_data(pd.DataFrame([ev.evaluate(inputs[0])])) G.generate(1)

Out[20]:

[{'x0': np.float64(-28.37156709762087), 'x1': np.float64(29.869929860151046)}]

In [21]:

# Create Xopt object
X = Xopt(evaluator=ev, vocs=vocs, generator=NelderMeadGenerator(vocs=vocs))

# Optional: give an initial pioint
X.generator.initial_point = {"x0": 0, "x1": 0}

# Create Xopt object X = Xopt(evaluator=ev, vocs=vocs, generator=NelderMeadGenerator(vocs=vocs)) # Optional: give an initial pioint X.generator.initial_point = {"x0": 0, "x1": 0}

In [22]:

X.run()

X.run()

In [23]:

# Generator is done and cannot be resumed
X.generator.is_done

# Generator is done and cannot be resumed X.generator.is_done

Out[23]:

True

In [24]:

# Generate calls will just return nothing
X.generator.generate(1) is None

# Generate calls will just return nothing X.generator.generate(1) is None

Out[24]:

True

In [25]:

# This shows the latest simplex
X.generator.simplex

# This shows the latest simplex X.generator.simplex

Out[25]:

{'x0': array([32.99996111, 32.99996171, 33.00002688]),
 'x1': array([44.00000851, 44.00006811, 44.00003045])}

In [26]:

X.data["y"].plot()
plt.yscale("log")

X.data["y"].plot() plt.yscale("log")

In [27]:

fig, ax = plt.subplots()
X.data.plot("x0", "x1", ax=ax, color="black", alpha=0.5)
ax.scatter(Xbest[0], Xbest[1], marker="x", color="red")

fig, ax = plt.subplots() X.data.plot("x0", "x1", ax=ax, color="black", alpha=0.5) ax.scatter(Xbest[0], Xbest[1], marker="x", color="red")

Out[27]:

<matplotlib.collections.PathCollection at 0x7f2f1131a710>

In [28]:

# This is the raw internal state of the generator
a = X.generator.current_state
a

# This is the raw internal state of the generator a = X.generator.current_state a

Out[28]:

SimplexState(astg=4, N=2, kend=3, jend=0, ind=array([2., 0., 1.]), sim=array([[32.99996111, 44.00000851],
       [32.99996171, 44.00006811],
       [33.00002688, 44.00003045]]), fsim=array([1.58463795e-09, 6.10453387e-09, 1.64942323e-09]), fxr=3.1653753661675e-08, x=array(nan), xr=array([32.99983049, 44.00005403]), xe=array([30.09830487, 46.28401515]), xc=array([32.99992966, 44.00013909]), xcc=array([33.00002688, 44.00003045]), xbar=array([32.99996141, 44.00003831]), doshrink=0, ngen=176)

In [29]:

# Check JSON representation of options
X.generator.json()

# Check JSON representation of options X.generator.json()

Out[29]:

'{"initial_point":{"x0":0.0,"x1":0.0},"initial_simplex":null,"adaptive":true,"xatol":0.0001,"fatol":0.0001,"current_state":{"astg":4,"N":2,"kend":3,"jend":0,"ind":[2.0,0.0,1.0],"sim":[[32.99996111240644,44.000008508408456],[32.99996171260399,44.0000681073357],[33.00002687564558,44.00003044869291]],"fsim":[1.5846379474430505e-9,6.104533869326014e-9,1.6494232253860804e-9],"fxr":3.1653753661675e-8,"x":null,"xr":[32.99983048622448,44.0000540262304],"xe":[30.098304873214147,46.28401515040761],"xc":[32.99992965625605,44.00013908571843],"xcc":[33.00002687564558,44.00003044869291],"xbar":[32.999961412505215,44.000038307872074],"doshrink":0,"ngen":176},"future_state":null,"x":null,"y":null,"is_done_bool":true}'

In [30]:

# Set the initial simplex to be the latest
X2 = Xopt(
    evaluator=ev,
    vocs=vocs,
    generator=NelderMeadGenerator(vocs=vocs, initial_simplex=X.generator.simplex),
)
X2.generator.xatol = 1e-9
X2.generator.fatol = 1e-9
X2.run()

X2.data["y"].plot()
plt.yscale("log")

# Set the initial simplex to be the latest X2 = Xopt( evaluator=ev, vocs=vocs, generator=NelderMeadGenerator(vocs=vocs, initial_simplex=X.generator.simplex), ) X2.generator.xatol = 1e-9 X2.generator.fatol = 1e-9 X2.run() X2.data["y"].plot() plt.yscale("log")

5-dimensional Rosenbrock¶

evaluate_rosenbrock works for arbitrary dimensions, so adding more variables to vocs transforms this problem.

In [31]:

YAML = """
generator:
  name: neldermead
evaluator:
  function: xopt.resources.test_functions.rosenbrock.evaluate_rosenbrock
vocs:
  variables:
    x1: [-5, 5]
    x2: [-5, 5]
    x3: [-5, 5]
    x4: [-5, 5]
    x5: [-5, 5]
  objectives:
    y: MINIMIZE
"""
X = Xopt.from_yaml(YAML)

YAML = """ generator: name: neldermead evaluator: function: xopt.resources.test_functions.rosenbrock.evaluate_rosenbrock vocs: variables: x1: [-5, 5] x2: [-5, 5] x3: [-5, 5] x4: [-5, 5] x5: [-5, 5] objectives: y: MINIMIZE """ X = Xopt.from_yaml(YAML)

In [32]:

X.run()
X.data["y"].plot()
plt.yscale("log")

X.run() X.data["y"].plot() plt.yscale("log")

In [33]:

fig, ax = plt.subplots(figsize=(8, 8))

Xgrid, Ygrid = np.meshgrid(np.linspace(-2, 2, 201), np.linspace(-2, 2, 201))

Zgrid = np.vectorize(lambda x, y: rosenbrock([x, y, 1, 1, 1]))(
    Xgrid, Ygrid
)  # The minimum is at 1,1,1,1,1
Zgrid = np.log(Zgrid + 1)

ax.pcolormesh(Xgrid, Ygrid, Zgrid)
ax.contour(Xgrid, Ygrid, Zgrid, levels=10, colors="black")
ax.set_xlabel("x0")
ax.set_ylabel("x1")


# Add all evaluations
ax.plot(X.data["x1"], X.data["x2"], color="red", alpha=0.5, marker=".")
ax.scatter(XMIN[0], XMIN[1], 50, marker="o", color="orange", label="True minimum")
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
# plt.legend()
ax.set_title("Xopt's Nelder-Mead progression")

fig, ax = plt.subplots(figsize=(8, 8)) Xgrid, Ygrid = np.meshgrid(np.linspace(-2, 2, 201), np.linspace(-2, 2, 201)) Zgrid = np.vectorize(lambda x, y: rosenbrock([x, y, 1, 1, 1]))( Xgrid, Ygrid ) # The minimum is at 1,1,1,1,1 Zgrid = np.log(Zgrid + 1) ax.pcolormesh(Xgrid, Ygrid, Zgrid) ax.contour(Xgrid, Ygrid, Zgrid, levels=10, colors="black") ax.set_xlabel("x0") ax.set_ylabel("x1") # Add all evaluations ax.plot(X.data["x1"], X.data["x2"], color="red", alpha=0.5, marker=".") ax.scatter(XMIN[0], XMIN[1], 50, marker="o", color="orange", label="True minimum") ax.set_xlim(-2, 2) ax.set_ylim(-2, 2) # plt.legend() ax.set_title("Xopt's Nelder-Mead progression")

Out[33]:

Text(0.5, 1.0, "Xopt's Nelder-Mead progression")