This is one of the 100 recipes of the IPython Cookbook, the definitive guide to high-performance scientific computing and data science in Python.
import numpy as np
import scipy as sp
import scipy.optimize as opt
import matplotlib.pyplot as plt
%matplotlib inline
f = lambda _: 1-np.sin(_)/_
x = np.linspace(-20., 20., 1000)
y = f(x)
plt.figure(figsize=(5,5));
plt.plot(x, y);
scipy.optimize
module comes with many function minimization routines. The minimize
function offers a unified interface to many algorithms. The Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm (default in minimize
) gives good results in general. The minimize
function requires an initial point as argument. For scalar univariate functions, we can also use minimize_scalar
.x0 = 3
xmin = opt.minimize(f, x0).x
Starting from $x_0=3$, the algorithm was able to find the actual global minimum, as shown on the following figure.
plt.figure(figsize=(5,5));
plt.plot(x, y);
plt.scatter(x0, f(x0), marker='o', s=300);
plt.scatter(xmin, f(xmin), marker='v', s=300);
plt.xlim(-20, 20);
x0 = 10
xmin = opt.minimize(f, x0).x
plt.figure(figsize=(5,5));
plt.plot(x, y);
plt.scatter(x0, f(x0), marker='o', s=300);
plt.scatter(xmin, f(xmin), marker='v', s=300);
plt.xlim(-20, 20);
xmin = opt.minimize(f, x0, method='Anneal').x
plt.figure(figsize=(5,5));
plt.plot(x, y);
plt.scatter(x0, f(x0), marker='o', s=300);
plt.scatter(xmin, f(xmin), marker='v', s=300);
plt.xlim(-20, 20);
This time, the algorithm was able to find the global minimum.
This function is very irregular and may be difficult to minimize in general. It is one of the many test functions for optimization that researchers have developed to study and benchmark optimization algorithms. (http://en.wikipedia.org/wiki/Test_functions_for_optimization)
def g(X):
# X is a 2*N matrix, each column contains
# x and y coordinates.
x, y = X
return np.sin(3*np.pi*x)**2+(x-1)**2*(1+np.sin(3*np.pi*y)**2)+(y-1)**2*(1+np.sin(2*np.pi*y)**2)
imshow
, on the square $[-10,10]^2$.n = 200
k = 10
X, Y = np.mgrid[-k:k:n*1j,-k:k:n*1j]
Z = g(np.vstack((X.ravel(), Y.ravel()))).reshape(n,n)
plt.figure(figsize=(5, 5));
# We use a logarithmic scale for the color here.
plt.imshow(np.log(Z), cmap=plt.cm.hot_r);
plt.xticks([]); plt.yticks([]);
x0, y0 = opt.minimize(g, (8, 3)).x
plt.figure(figsize=(5, 5));
plt.imshow(np.log(Z), cmap=plt.cm.hot_r,
extent=(-k, k, -k, k), origin=0);
plt.scatter(x0, y0, s=100);
plt.xticks([]); plt.yticks([]);
You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).
IPython Cookbook, by Cyrille Rossant, Packt Publishing, 2014 (500 pages).