Rosenbrock banana

Here, we perform optimization for the Rosenbrock banana function, which does not require an AMICI model. In particular, we try several ways of specifying derivative information.

[1]:
import pypesto
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

%matplotlib inline

Define the objective and problem

[2]:
# first type of objective
# get the optimizer trace
objective_options = pypesto.ObjectiveOptions(trace_record=True, trace_save_iter=1)
objective1 = pypesto.Objective(fun=sp.optimize.rosen,
                               grad=sp.optimize.rosen_der,
                               hess=sp.optimize.rosen_hess,
                               options=objective_options)

# second type of objective
def rosen2(x):
    return sp.optimize.rosen(x), sp.optimize.rosen_der(x), sp.optimize.rosen_hess(x)
objective2 = pypesto.Objective(fun=rosen2, grad=True, hess=True)

dim_full = 10
lb = -5 * np.ones((dim_full, 1))
ub = 5 * np.ones((dim_full, 1))

problem1 = pypesto.Problem(objective=objective1, lb=lb, ub=ub)
problem2 = pypesto.Problem(objective=objective2, lb=lb, ub=ub)

Illustration

[3]:
x = np.arange(-2, 2, 0.1)
y = np.arange(-2, 2, 0.1)
x, y = np.meshgrid(x, y)
z = np.zeros_like(x)
for j in range(0, x.shape[0]):
    for k in range(0, x.shape[1]):
        z[j,k] = objective1([x[j,k], y[j,k]], (0,))
[4]:
fig = plt.figure()
fig.set_size_inches(*(14,10))
ax = plt.axes(projection='3d')
ax.plot_surface(X=x, Y=y, Z=z)
plt.xlabel('x axis')
plt.ylabel('y axis')
ax.set_title('cost function values')
[4]:
Text(0.5, 0.92, 'cost function values')
../_images/example_rosenbrock_7_1.png

Run optimization

[5]:
# create different optimizers
optimizer_bfgs = pypesto.ScipyOptimizer(method='l-bfgs-b')
optimizer_tnc = pypesto.ScipyOptimizer(method='TNC')
optimizer_dogleg = pypesto.ScipyOptimizer(method='dogleg')

# set number of starts
n_starts = 20

# Run optimizaitons for different optimzers
result1_bfgs = pypesto.minimize(problem=problem1, optimizer=optimizer_bfgs, n_starts=n_starts)
result1_tnc = pypesto.minimize(problem=problem1, optimizer=optimizer_tnc, n_starts=n_starts)
result1_dogleg = pypesto.minimize(problem=problem1, optimizer=optimizer_dogleg, n_starts=n_starts)

# Optimize second type of objective
result2 = pypesto.minimize(problem=problem2, optimizer=optimizer_tnc, n_starts=n_starts)
/home/yannik/yenv/lib/python3.6/site-packages/scipy/optimize/_minimize.py:518: RuntimeWarning: Method dogleg cannot handle constraints nor bounds.
  RuntimeWarning)

Visualize and compare optimization results

[6]:
import pypesto.visualize

# plot separated waterfalls
pypesto.visualize.waterfall(result1_bfgs, size=(15,6))
pypesto.visualize.waterfall(result1_tnc, size=(15,6))
pypesto.visualize.waterfall(result1_dogleg, size=(15,6))
[6]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fac6554b550>
../_images/example_rosenbrock_11_1.png
../_images/example_rosenbrock_11_2.png
../_images/example_rosenbrock_11_3.png

We can now have a closer look, which method perfomred better: Let’s first compare bfgs and TNC, since both methods gave good results. How does they fine convergence look like?

[7]:
# plot one list of waterfalls
pypesto.visualize.waterfall([result1_bfgs, result1_tnc],
                            legends=['L-BFGS-B', 'TNC'],
                            start_indices=10,
                            scale_y='lin')
[7]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fac65342438>
../_images/example_rosenbrock_13_1.png
[8]:
# retrieve second optimum
all_x = result1_bfgs.optimize_result.get_for_key('x')
all_fval = result1_bfgs.optimize_result.get_for_key('fval')
x = all_x[19]
fval = all_fval[19]
print('Second optimum at: ' + str(fval))

# create a reference point from it
ref = {'x': x, 'fval': fval, 'color': [
    0.2, 0.4, 1., 1.], 'legend': 'second optimum'}
ref = pypesto.visualize.create_references(ref)

# new waterfall plot with reference point for second optimum
pypesto.visualize.waterfall(result1_dogleg, size=(15,6),
                            scale_y='lin', y_limits=[-1, 101],
                            reference=ref, colors=[0., 0., 0., 1.])
Second optimum at: 3.9865791128196464
[8]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fac652265f8>
../_images/example_rosenbrock_14_2.png

Visualize parameters

There seems to be a second local optimum. We want to see whether it was also found by the dogleg method

[9]:
pypesto.visualize.parameters([result1_bfgs, result1_tnc],
                            legends=['L-BFGS-B', 'TNC'],
                            balance_alpha=False)
pypesto.visualize.parameters(result1_dogleg,
                             legends='dogleg',
                             reference=ref,
                             size=(15,10),
                             start_indices=[0, 1, 2, 3, 4, 5],
                             balance_alpha=False)
[9]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fac654570b8>
../_images/example_rosenbrock_17_1.png
../_images/example_rosenbrock_17_2.png

Optimizer history

Let’s compare optimzer progress over time.

[10]:
# plot one list of waterfalls
pypesto.visualize.optimizer_history([result1_bfgs, result1_tnc],
                                    legends=['L-BFGS-B', 'TNC'],
                                    reference=ref)
# plot one list of waterfalls
pypesto.visualize.optimizer_history(result1_dogleg,
                                    reference=ref)
[10]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fac6549ef60>
../_images/example_rosenbrock_20_1.png
../_images/example_rosenbrock_20_2.png

We can also visualize this usign other scalings or offsets…

[11]:
# plot one list of waterfalls
pypesto.visualize.optimizer_history([result1_bfgs, result1_tnc],
                                    legends=['L-BFGS-B', 'TNC'],
                                    reference=ref,
                                    offset_y=0.)

# plot one list of waterfalls
pypesto.visualize.optimizer_history([result1_bfgs, result1_tnc],
                                    legends=['L-BFGS-B', 'TNC'],
                                    reference=ref,
                                    scale_y='lin',
                                    y_limits=[-1., 11.])
[11]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fac6501beb8>
../_images/example_rosenbrock_22_1.png
../_images/example_rosenbrock_22_2.png

Compute profiles

The profiling routine needs a problem, a results object and an optimizer.

Moreover it accepts an index of integer (profile_index), whether or not a profile should be computed.

Finally, an integer (result_index) can be passed, in order to specify the local optimum, from which profiling should be started.

[12]:

# compute profiles
profile_options = pypesto.ProfileOptions(min_step_size=0.0005,
    delta_ratio_max=0.05,
    default_step_size=0.005,
    ratio_min=0.03)

result1_tnc = pypesto.parameter_profile(
    problem=problem1,
    result=result1_tnc,
    optimizer=optimizer_tnc,
    profile_index=np.array([1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0]),
    result_index=0,
    profile_options=profile_options)

# compute profiles from second optimum
result1_tnc = pypesto.parameter_profile(
    problem=problem1,
    result=result1_tnc,
    optimizer=optimizer_tnc,
    profile_index=np.array([1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0]),
    result_index=19,
    profile_options=profile_options)

Visualize and analyze results

pypesto offers easy-to-use visualization routines:

[13]:
# specify the parameters, for which profiles should be computed
ax = pypesto.visualize.profiles(result1_tnc, profile_indices = [0,1,2,5],
                           reference=ref, profile_list_id=0)
# plot profiles again, now from second optimum
ax = pypesto.visualize.profiles(result1_tnc, profile_indices = [0,1,2,5],
                           reference=ref, profile_list_id=1)
../_images/example_rosenbrock_28_0.png
../_images/example_rosenbrock_28_1.png

If the result needs to be examined in more detail, it can easily be exported as a pandas.DataFrame:

[14]:
result1_tnc.optimize_result.as_dataframe(['fval', 'n_fval', 'n_grad',
                                        'n_hess', 'n_res', 'n_sres', 'time'])
[14]:
fval n_fval n_grad n_hess n_res n_sres time
0 3.727426e-13 204 204 0 0 0 2.277306
1 4.994350e-13 200 200 0 0 0 2.517795
2 2.106192e-12 207 207 0 0 0 2.310546
3 2.368548e-12 226 226 0 0 0 2.805628
4 2.407578e-12 95 95 0 0 0 1.029478
5 5.525261e-12 166 166 0 0 0 1.876049
6 1.038073e-11 249 249 0 0 0 2.957421
7 3.148873e-11 216 216 0 0 0 2.358464
8 9.991588e-11 170 170 0 0 0 2.292490
9 2.045469e-09 156 156 0 0 0 2.203442
10 2.414312e-09 187 187 0 0 0 2.385922
11 3.400104e-09 264 264 0 0 0 3.420013
12 4.277078e-09 158 158 0 0 0 2.009642
13 6.561378e-09 171 171 0 0 0 1.820322
14 2.038944e-08 181 181 0 0 0 2.000105
15 3.986579e+00 158 158 0 0 0 1.917194
16 3.986579e+00 106 106 0 0 0 1.417423
17 3.986579e+00 90 90 0 0 0 1.017045
18 3.986579e+00 138 138 0 0 0 1.473538
19 3.986579e+00 155 155 0 0 0 2.147978