# SciPy_optimization.ipynb

notebooks/SciPy_optimization.ipynb

# Notebook to demonstrate SciPY optimization methods¶

### Dr. Tirthajyoti Sarkar, Fremont, CA 94536¶

Mathematical optimization is at the heart of solutions to major business problems in engineering, finance, healthcare, socioeconomic affairs. Pretty much all business problems boil down to minimization of some kind of resource cost or maximization of some kind of profit given other constraints.

An optimization process is also the soul of operation research, which is intimately related to modern data-driven business analytics. In this manner, it is also closely related to the data science pipeline, employed in virtually all businesses today.Â

Although much has been written about the data wrangling and predictive modeling aspects of a data science project, the final frontier often involves solving an optimization problem using the data-driven models which can improve the bottom-line of the business by reducing cost or enhancing productivity.

Python has become the de-facto lingua franca of analytics, data science, and machine learning. Therefore, it makes sense to discuss optimization packages and frameworks within the Python ecosystem.

We cover optimization algorithms available within the SciPy ecosystem. SciPy is the most widely used Python package for scientific and mathematical analysis and it is no wonder that it boasts of powerful yet easy-to-use optimization routines for solving complex problems.

#### SciPy optimization reference guide¶

import numpy as np
from matplotlib import pyplot as plt
from scipy import optimize

### Minimize a simple scalar function: $\text{sin}(x).\text{exp}[(x-0.6)^2]$¶

def scalar1(x):
return np.sin(x)*np.exp(-0.1*(x-0.6)**2)
def plot_nice(x,y,title=None,xlabel='x',ylabel='y',show=True):
#plt.figure(figsize=(8,5))
if title!=None:
plt.title(str(title)+'\n',fontsize=18)
plt.plot(x,y,color='k',lw=3)
plt.grid(True)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.xlabel(xlabel,fontsize=15)
plt.ylabel(ylabel,fontsize=15)
if show:
plt.show()
x = np.arange(-10,10,0.05)
y = scalar1(x)
plot_nice(x,y,title="Objective function",xlabel='x-values',ylabel='Function values')

### Use optimize.minimize_scalar() method¶

result = optimize.minimize_scalar(scalar1)
result['success']
True
print("Minimum occurs at: ",result['x'])
Minimum occurs at:  -1.2214484245210282
print(result)
fun: -0.6743051024666711
nfev: 15
nit: 10
success: True
x: -1.2214484245210282

### Bounded search (bound on the independent variable)¶

result = optimize.minimize_scalar(scalar1,bounds=(0,10),method='Bounded')
print("When bounded between 0 and 10, minimum occurs at: ",result['x'])
When bounded between 0 and 10, minimum occurs at:  4.101466164987216

### Other function-based constraints¶

def constraint1(x):
return 0.5-np.log10(x**2+2)
def constraint2(x):
return np.log10(x**2+2) - 1.5
def constraint3(x):
return np.sin(x)+0.3*x**2-1
cons = ({'type':'ineq','fun':constraint1},
{'type':'ineq','fun':constraint2},
{'type':'eq','fun':constraint3})
result = optimize.minimize(scalar1,x0=0,method='SLSQP',constraints=cons,options={'maxiter':100})
print(result)
fun: 0.7631695862891654
jac: array([0.59193639])
message: 'Iteration limit exceeded'
nfev: 1241
nit: 101
njev: 100
status: 9
success: False
x: array([0.8773752])
result = optimize.minimize(scalar1,x0=-20,method='SLSQP',constraints=cons,options={'maxiter':100})
print(result)
fun: -0.28594944567686104
jac: array([-0.46750661])
message: 'Iteration limit exceeded'
nfev: 1233
nit: 101
njev: 101
status: 9
success: False
x: array([-2.37569791])
result = optimize.minimize(scalar1,x0=-20,method='SLSQP',constraints=cons,options={'maxiter':3})
print(result)
fun: -0.4155114388552631
jac: array([-0.46860977])
message: 'Iteration limit exceeded'
nfev: 12
nit: 4
njev: 4
status: 9
success: False
x: array([-2.10190632])

### Multi-variate case: Sum of Gaussians¶

mu = [-1,0.3,2.1]
sigma = [2.1,0.8,1.7]
plt.figure(figsize=(8,5))
for m,s in zip(mu,sigma):
x=np.arange(-5,5,0.01)
y=np.exp(-(x-m)**2/(2*s**2))
plot_nice(x,y,show=False)

plt.show()
def gaussian(m,s):
x = np.arange(-5,5,0.01)
return np.exp(-(x-m)**2/(2*s**2))
def gaussian_mixture(x):
"""
Computes the resultant Gaussian mixture from an input vector and known mean, variance quantities
"""
return -(np.exp(-(x[0]+1)**2/(2.1**2))+np.exp(-(x[1]-0.3)**2/(0.8**2))+np.exp(-(x[2]-2.1)**2/(1.7**2)))
gaussian_mixture(np.array([3,-2,2]))
-1.0233691172715331
x0=np.array([0]*3)
result = optimize.minimize(gaussian_mixture,x0=x0,method='SLSQP',options={'maxiter':100})
result
fun: -2.9999996182639146
jac: array([-8.09729099e-05, -2.40176916e-04,  7.11053610e-04])
message: 'Optimization terminated successfully.'
nfev: 42
nit: 8
njev: 8
status: 0
success: True
x: array([-1.00017856,  0.29992313,  2.10102744])

### Bounds with multiple variables¶

x0=np.array([0]*3)
x1_bound = (-2,2)
x2_bound = (0,5)
x3_bound = (-3,0)
result = optimize.minimize(gaussian_mixture,x0=x0,method='SLSQP',options={'maxiter':100},
bounds=(x1_bound,x2_bound,x3_bound))
result
fun: -2.217414055755018
jac: array([-2.89082527e-06,  3.60012054e-04, -3.15965086e-01])
message: 'Optimization terminated successfully.'
nfev: 31
nit: 6
njev: 6
status: 0
success: True
x: array([-1.00000644e+00,  3.00115191e-01, -8.03574200e-17])