Problem 4¶
What is the global minimum of the function $ f(x,y) = e^{\sin(50x)} + \sin(60e^y) + \sin(70\sin{x}) + \sin(\sin{80y}) - \sin(10(x+y)) + (x^2+y^2)/4 $ ?
Plotting the function given, we note that $f$ is bounded by the $(x^2 + y^2)/4$ term for large $x,y$; however, the function has many local minima/maxima at $x,y \approx 0$ making it difficult to find a global minimum.
#from math import exp, sin
import numpy as np
from numpy import exp, sin
import matplotlib.pyplot as plt
from matplotlib import cm
f = lambda x,y: exp(sin(50*x)) + sin(60*exp(y)) + sin(70*sin(x)) + sin(sin(80*y)) - sin(10*(x+y)) + (x*x+y*y)/4
X = np.arange(-10,10,0.5)
Y = np.arange(-10,10,0.5)
X, Y = np.meshgrid(X,Y)
Z = f(X,Y)
# fig = plt.figure()
fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.view_init(elev=10)
surf = ax.plot_surface(X,Y,Z,cmap=cm.hsv,linewidth=20)
X = np.arange(-0.5,0.5,0.001)
Y = np.arange(-0.5,0.5,0.001)
X, Y = np.meshgrid(X,Y)
Z = f(X,Y)
ax = fig.add_subplot(1, 2, 2, projection='3d')
surf = ax.plot_surface(X,Y,Z,cmap=cm.coolwarm,linewidth=20,antialiased=False)
Paper and pencil minimization¶
Note that we have
$ \begin{cases} \min(e^{\sin(50x)}) = 1/e\\ \min(\sin(60 e^y)) = -1 \\ \min(\sin(\sin 70x)) = -\sin1 \\ \min(\sin(\sin 80y)) = -\sin1 \\ \min(-\sin(10(x+y)) = -1 \\ \min((x^2+y^2)/4) = 0 \\ \end{cases} \implies \min(f(x,y)) \approx 1/e-2-2\sin1 \approx -3.315 $
Using built-in methods¶
f = lambda x,y: exp(sin(50*x)) + sin(60*exp(y)) + sin(70*sin(x)) + sin(sin(80*y)) - sin(10*(x+y)) + (x*x+y*y)/4
f_min = lambda x: f(x[0],x[1]) # wrapper for f so that it can be passed to scipy.optimize
from scipy.optimize import minimize
from time import time
methods = ["Nelder-Mead","Powell","CG","BFGS","L-BFGS-B","TNC","COBYLA","SLSQP","trust-constr"]
def time_method(method): # wrapper function to help time scipy.optimize with given method
start = time()
out = minimize(f_min,(0,0),method=method,tol=1e-12)
return (out.fun, out.x, time()-start, method)
for method in methods:
print(time_method(method))
(-0.7130747093105191, array([-0.02230223, -0.00472762]), 0.007046937942504883, 'Nelder-Mead') (-2.966665564173802, array([ 0.24280684, -0.09332387]), 0.011039495468139648, 'Powell') (-2.0990843569596604, array([-1.01123843, -0.09335979]), 0.01640629768371582, 'CG') (-2.0990843569595934, array([-1.01123844, -0.09335979]), 0.01039576530456543, 'BFGS') (-2.0990843569597404, array([-1.01123844, -0.09335979]), 0.0033600330352783203, 'L-BFGS-B') (-0.7130747093103592, array([-0.02230224, -0.00472763]), 0.006999015808105469, 'TNC') (-0.8627887931635007, array([-1.65508048, 0.51042469]), 0.0015561580657958984, 'COBYLA') (-0.6254386181940512, array([-0.18837361, 0.77441795]), 0.004513263702392578, 'SLSQP') (-1.452607474505165, array([-1.16099618, -0.00550989]), 0.06516742706298828, 'trust-constr')
As we can see, none of these methods yield well to our given function. The Powell method provides the lowest value of the functio, however, calculating the derivative of the function at the point given:
from scipy.optimize import approx_fprime
print(approx_fprime([ 0.24280684, -0.09332387], f_min, [1e-12,1e-12]))
[-0.00177635 0. ]
we have that that the partial derivative with respect to $x$ of the function at this point is negative. Thus, we cannot have this as a global minimum of this function as we can show
print(f(0.2428064,-0.09332387)) # Minimum provided by scipy.optimize
print(f(0.2428066,-0.09332387)) # Lower value of fuction
-2.9666655635549026 -2.966665563987929
Bringing out the big guns¶
Since this function has many local minima, we have to resort to using more specialized algorithms in order to find a global minimum. We use two below
Basin-hopping¶
The basin-hopping algorithm described in a 1997 paper by David J Wales and Jonathan Doye is a global optimization algorithm originally developed for finding the minimum energy structure of molecules.
From the scipy documentation of the implementation of this function:
Designed to mimic the natural process of energy minimization of clusters of atoms, it works well for similar problems with “funnel-like, but rugged” energy landscapes
Thus, it should work well for our function:
from scipy.optimize import basinhopping
import time
start = time.time()
print(basinhopping(f_min, [0,0], niter = 5000))
print("Time elapsed: ", time.time()-start, " seconds")
fun: -3.3068686474749844
lowest_optimization_result: fun: -3.3068686474749844
hess_inv: array([[ 1.66496802e-04, -1.63442544e-06],
[-1.63442544e-06, 1.01069713e-04]])
jac: array([1.18911266e-05, 1.63316727e-05])
message: 'Desired error not necessarily achieved due to precision loss.'
nfev: 168
nit: 8
njev: 53
status: 2
success: False
x: array([-0.02440309, 0.21061242])
message: ['requested number of basinhopping iterations completed successfully']
minimization_failures: 980
nfev: 325678
nit: 5000
njev: 105109
success: False
x: array([-0.02440309, 0.21061242])
Time elapsed: 33.00804781913757 seconds
Differential Evolution¶
Differential evolution is a stochastic genetic algorithm that is often used for global optimization problems. Briefly, some initial candidates for the minimum of the function are generated. Then these candidates are "bred" and similarly to the process of natural selection, over multiple generations, candidates that are closer to the minimum of the function are generated
from scipy.optimize import differential_evolution
import time
start = time.time()
print(differential_evolution(f_min, [[-1,1],[-1,1]],tol=1e-12))
print("Time elapsed: ", time.time()-start, " seconds")
fun: -3.306868647475183
message: 'Optimization terminated successfully.'
nfev: 1956
nit: 63
success: True
x: array([-0.02440308, 0.21061243])
Time elapsed: 0.25775575637817383 seconds
Wrapping it up: more digits¶
Since we now have the location of the global minimum of the function, we can use local optimization algorithms to get more digits of the function's minimum. In this case, we'll start with the result given above as the minimum of the function and then do a simple grid search in the neighborhood of this point, updating the resulting minimum until we get the desired accuracy.
from mpmath import mp, mpf, exp, sin
from tqdm import tqdm
mp.dps = 1000
f = lambda x,y: exp(sin(50*x)) + sin(60*exp(y)) + sin(70*sin(x)) + sin(sin(80*y)) - sin(10*(x+y)) + (x*x+y*y)/4
x = mpf("-0.02440308")
y = mpf("0.21061242")
mins = [None]
step = mpf("1e-7")
for i in tqdm(range(10000)):
neighborhood = [[x,y+step,f(x,y+step)],
[x,y-step,f(x,y-step)],
[x+step,y,f(x+step,y)],
[x-step,y,f(x-step,y)],
[x-step,y-step,f(x-step,y-step)],
[x-step,y+step,f(x-step,y+step)],
[x+step,y-step,f(x+step,y-step)],
[x+step,y+step,f(x+step,y+step)]]
min_f_value = min(neighborhood, key = lambda x: x[2])
x,y,minimum = min_f_value
mins.append(minimum)
if mins[-1]==mins[-2]: print(i); break
step /= 2
print(mins[-1])
16%|██████████████████▉ | 1645/10000 [01:10<05:58, 23.30it/s]
1645 -3.306868647475237280076113770898515657166482361476288217501293085496309199837888295035825488075283499186192630891136430717132695479183570527535481485323493119666055496491770987805562739806177408729014170264272392227086207334639477551474403455597059000710130930484549790041257397877006760662403161597724611217893259536294053937533901074129148867476633116962279417362942396292658321079163381013168510632678554019900401524711714355688234420637521685431767309457909425768882351097148504853083520986326339983557947849154630104547099679467311070004309817910593652674429408839748633882190364165440517209476651683601753535141475760124791643629970356353417693895735350598635473642122822528688878749695509827898143628521390494115341245395385269727851257267362553188546604443742861873679538628389385868997963047779486318537165899052250504507672343475949688780439524509796737104261919981577990979508967297886317173392807994861717751825429547302477283761015961319823829914040695514974564405118425411424061122303324