%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy
import scipy.optimize as opt
import scipy.special as spec
import seaborn as sns
sns.set(
style='ticks',
context='talk'
)
Scipy has several methods for optimizing (minimizing or maximizing) functions.
This session is roughly based on Tutorial on "Modern Optimization Methods in Python" .
Consider the following polynomial, defined using np.poly1d1
, which, given the polynomial coeffcients $a_i$ returns a polynomial $\sum_{i=0}^n{a_i x^i}$:
p = np.poly1d((1.3, 4, 0.6))
print(p)
2 1.3 x + 4 x + 0.6
def minplot(x, f, xmin=None):
plt.plot(x, f(x))
if xmin is not None:
plt.plot(xmin, f(xmin), 'ok')
plt.axvline(xmin, color='k', ls='--')
plt.xlabel('x')
plt.ylabel('f(x)')
sns.despine()
x = np.linspace(-10, 7)
minplot(x, p)
We'll use scipy.optimize.fmin
(aliasing scipy.optimize
as opt
) to find a minimum for p
, our polynomial. fmin
uses the downhill simplex algorithm and we start with the guess 3:
xmin = opt.fmin(p, [3])
print("x = {0}, p(x) = {1}".format(xmin, p(xmin)))
minplot(x, p, xmin)
Optimization terminated successfully. Current function value: -2.476923 Iterations: 20 Function evaluations: 40 x = [-1.53845215], p(x) = [-2.47692308]
Find the minimum of the function $f(x) = (x-2)(x+1)^3$ using either fmin
or minimize_scalar
. The later doesn't require a guess (x0
).
f(1.250) = -8.543
Suppose we want to find the minimum of a 1st order Bessel function with the constraint $2<x<4$:
x = np.linspace(0, 20, 500)
minplot(x, spec.j1)
plt.axvspan(8, 16, alpha=0.3);
We use opt.minimize_scalar
and give it constraints:
result = opt.minimize_scalar(spec.j1, method="bounded", bounds=[8, 16])
print(result)
minplot(x, spec.j1, result.x)
plt.axvspan(8, 16, alpha=0.3);
fun: -0.23330441717143405 message: 'Solution found.' nfev: 9 status: 0 success: True x: 11.706004881285123
Find the maximum of the Fresnel sin integral, which is the first value returned from the $\mathbb R \to \mathbb R^2$ function scipy.special.fresnel
. Constrain maximization to $x \in (2, 3)$. Plot the result and the maximum to make sure you got it right.
Follows example from SciPy Tutorial.
The opt.minimize
provides interface to other contrained methods. Let's use SLSQP (Sequential Least SQuares Programming) method to deal with the following problem:
def f(arr, sign=1.0):
x, y = arr
return (2 * x * y + 2 * x - x**2 - 2 * y**2) * sign
x = np.linspace(-10, 10)
y = np.linspace(-10, 10)
We load the mplot3d
module that will allow us to do some 3D plotting.
from mpl_toolkits.mplot3d import Axes3D
def minplot3d(x, y, f, xmin=None, ymin=None):
z = np.array([[f([xi, yi]) for yi in y] for xi in x])
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, rstride=4, cstride=4, cmap='viridis', alpha=0.5)
if xmin is not None and ymin is not None:
ax.scatter3D(xmin, ymin, f([xmin, ymin]), s=50, c='k')
minplot3d(x, y, f)
For this method we also need the derivative of f
with respect to x
and y
:
def df(arr, sign=1.0):
x, y = arr
dfdx = -2*x + 2*y + 2
dfdy = 2*x - 4*y
return np.array([ dfdx, dfdy ]) * sign
minplot3d(x, y, lambda x: df(x)[0])
minplot3d(x, y, lambda x: df(x)[1])
Now an unconstrained optimization can be performed as:
res = opt.minimize(
f,
x0=[-1, 1],
args=(-1.0,),
jac=df,
method='SLSQP',
options={'disp': True}
)
print(res)
print("Min at ", res.x)
xmin, ymin = res.x
minplot3d(x, y, f, xmin, ymin)
Optimization terminated successfully. (Exit mode 0) Current function value: -2.0 Iterations: 4 Function evaluations: 5 Gradient evaluations: 4 fun: -2.0 jac: array([-0., -0.]) message: 'Optimization terminated successfully.' nfev: 5 nit: 4 njev: 4 status: 0 success: True x: array([2., 1.]) Min at [2. 1.]
Then constraints are defined as a sequence of dictionaries, with keys type
, fun
and jac
. Reminder:
If we had a $<=$ inequality, we would have multiplied both sides by -1.
cons = (
dict(
type='eq',
fun=lambda v: np.array([v[0]**3 - v[1]]), # v[0] -> x, v[1] -> y
jac=lambda v: np.array([3 * (v[0]**2), -1])
),
dict(
type='ineq',
fun=lambda v: np.array([v[1] - 1]),
jac=lambda v: np.array([0, 1])
)
)
res = opt.minimize(
f,
x0=[-1, 1],
args=(-1.0,),
jac=df,
constraints=cons,
method='SLSQP',
options={'disp': True}
)
print(res)
print()
print("Min at ", res.x)
xmin, ymin = res.x
print(xmin**3 - ymin)
print(ymin - 1)
Optimization terminated successfully. (Exit mode 0) Current function value: -1.0000001831052137 Iterations: 9 Function evaluations: 14 Gradient evaluations: 9 fun: -1.0000001831052137 jac: array([-1.99999982, 1.99999982]) message: 'Optimization terminated successfully.' nfev: 14 nit: 9 njev: 9 status: 0 success: True x: array([1.00000009, 1. ]) Min at [1.00000009 1. ] 2.746578582346615e-07 0.0
minplot3d(x, y, f, xmin, ymin)
Minimize the Rosenbrock function:
$$
f(x) = \sum_{i=1}^{N-1}{100(x_i - x_{i-1}^2)^2 + (1-x_{i-1})^2}
$$
defined in opt.rosen
as a function that accepts an array x
.
Use how many dimensions you want (determined by the initial guess), use the function opt.minimize
with method
set to either nelder-mead
or powell
.
The result should be an array of 1
s. If you get something close but not quite, you can change the tolerance by giving opt.minimize
the argument options={'xtol':tol}
, where tol
is the error tolerance you want.
This notebook was written by Yoav Ram and is part of the Python for Engineers course.
The notebook was written using Python 3.6.1. Dependencies listed in environment.yml, full versions in environment_full.yml.
This work is licensed under a CC BY-NC-SA 4.0 International License.