07.01 - OPTIMIZATION

!wget --no-cache -O init.py -q https://raw.githubusercontent.com/fagonzalezo/ai4eng-unal/main/content/init.py
import init; init.init(force_download=False); init.get_weblink()

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import plotly.graph_objects as go
%matplotlib inline

def plot3D(X0, X1, Z, initx=None, x=None, hist=None):
    data=[ go.Surface(z=Z, x=X0, y=X1, opacity=.5) ]
    if initx is not None:
        data.append(go.Scatter3d(x=[initx[0]], y=[initx[1]], z=[f(*initx)], mode='markers',
                                         marker_symbol="x",
                                          marker=dict(size=4,opacity=1, color="black"), name="initx"))
        
    if x is not None:
        data.append(go.Scatter3d(x=[x[0]], y=[x[1]], z=[f(*x)], mode='markers', name="min",
                                          marker=dict(size=4,opacity=1, color="black")))        
        
        
    if hist is not None:
        data.append(go.Scatter3d(x=hist[:,0], y=hist[:,1], z=f(hist[:,0], hist[:,1]), 
                    marker_symbol="x", marker=dict(size=2,opacity=0.5, color="black"),
                    line=dict(color='darkblue',width=2), name="optimization"))
        
    fig = go.Figure(data=data)
    fig.update_layout(autosize=False, width=500, height=500)
    fig.update_layout(legend=dict(x=-.1, y=1))
    fig.show()  

Black box optimization

Given a function

\[f(x)=(x-3)^2+2\]

We want to find the minimum value of the function

\[\underset{x}{\text{min}} \;\;f(x)\]

and the value of \(x\) at which the minimum value is reached

\[\underset{x}{\text{arg min}} \;\;f(x)\]

Observe how we can use black box optimization algorithms from scipy.optimize.minimize.

from scipy.optimize import minimize

f = lambda x: (x-3)**2+2

r = minimize(f, method="BFGS", x0=np.random.normal()*10 ) 
r

      fun: 2.0
 hess_inv: array([[0.5]])
      jac: array([0.])
  message: 'Optimization terminated successfully.'
     nfev: 15
      nit: 4
     njev: 5
   status: 0
  success: True
        x: array([3.])

## KEEPOUTPUT
xr = np.linspace(-1,5)
plt.plot(xr, f(xr), color="red", lw=3)
plt.grid();
plt.axvline(r.x, color="black", label=r"$arg min \;\;f(x)=%.2f$"%r.x)
plt.axhline(r.fun, color="black", ls="--", label=r"$min \;\;f(x)=%.2f$"%r.fun)
plt.legend()

<matplotlib.legend.Legend at 0x7f5c2569ec70>

the minimum might be local

f = lambda x: (x-3)**2*np.sin(x)**2+x

initx = np.random.random()*6-1
r = minimize(f, method="BFGS", x0=initx) 
r

      fun: -0.02684262682443163
 hess_inv: array([[0.0503387]])
      jac: array([6.9453381e-07])
  message: 'Optimization terminated successfully.'
     nfev: 24
      nit: 6
     njev: 8
   status: 0
  success: True
        x: array([-0.05283396])

## KEEPOUTPUT
xr = np.linspace(np.min([-1,r.x-1]), np.max([5, r.x+1]))
plt.plot(xr, f(xr), color="red", lw=3)
plt.grid();
plt.axvline(r.x, color="black", label=r"$arg min \;\;f(x)=%.2f$"%r.x)
plt.axhline(r.fun, color="black", ls="--", label=r"$min \;\;f(x)=%.2f$"%r.fun)
plt.scatter(initx, f(initx), color="black", marker="x", lw=5, s=150, label="initx")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))

<matplotlib.legend.Legend at 0x7f5c259045e0>

minimizing a multivariable function

Given a function

\[f(x)=(x*y-3)^2+2\]

We want to find the minimum value of the function

\[\underset{x_0,x_1}{\text{min}} \;\;f(x_0, x_1)\]

and the value of \(x\) at which the minimum value is reached

\[\underset{x_0,x_1}{\text{arg min}} \;\;f(x_0, x_1)\]

f = lambda x0,x1: x0**2+x1**2

initx = np.random.random(size=2)*10-5
r = minimize(lambda x: f(*x), method="BFGS", x0=initx)
r

      fun: 1.1163460495926846e-16
 hess_inv: array([[ 0.69131361, -0.24301421],
       [-0.24301421,  0.80868639]])
      jac: array([ 2.79045959e-10, -3.54400608e-10])
  message: 'Optimization terminated successfully.'
     nfev: 16
      nit: 2
     njev: 4
   status: 0
  success: True
        x: array([-7.31105762e-09, -7.62778090e-09])

x0 = np.linspace(-5, 5, 100)
x1 = np.linspace(-5, 5, 100)

X0, X1 = np.meshgrid(x0, x1)
Z = f(X0, X1)

## KEEPOUTPUT
plot3D(X0,X1,Z, initx, r.x)

observe how we use vectorization out of the box, and a meshgrid with all possible input values for \(x_0\) and \(x_1\)

## KEEPOUTPUT
x0.shape, x1.shape, X0.shape, X1.shape, Z.shape

((100,), (100,), (100, 100), (100, 100), (100, 100))

## KEEPOUTPUT
plt.figure(figsize=(8,3))
plt.subplot(121); plt.contourf(X0, X1, X0, levels=len(x0), cmap=plt.cm.coolwarm); plt.colorbar(); plt.title("X0")
plt.subplot(122); plt.contourf(X0, X1, X1, levels=len(x1), cmap=plt.cm.coolwarm); plt.colorbar(); plt.title("X1")

Text(0.5, 1.0, 'X1')

f = lambda x0,x1: np.cos(.5*x1)+np.sin(.5*x0)

sigm = lambda x: 1/(1+np.exp(-x))
f = lambda x0,x1: sigm(2*np.sin(2*x0)*np.cos(2*x1)*(.1*x0+1)**2)

def log(xk):
    h.append(xk)
    
initx = np.random.random(size=2)*10-5
h = [initx]
r = minimize(lambda x: f(*x), method="BFGS", x0=initx, callback=log)
h = np.r_[h]
r

      fun: 0.1534065009294887
 hess_inv: array([[ 1.05619925, -0.02483636],
       [-0.02483636,  1.11495152]])
      jac: array([-1.53295696e-06, -7.30156898e-07])
  message: 'Optimization terminated successfully.'
     nfev: 40
      nit: 8
     njev: 10
   status: 0
  success: True
        x: array([-7.31660667e-01, -8.27436002e-07])

np.min(h, axis=0), np.max(h, axis=0)

(array([-0.73194762, -1.02659659]), array([-0.41815349,  0.0229568 ]))

## KEEPOUTPUT
x0 = np.linspace(np.min(h[:,0]-1), np.max(h[:,0]+1), 50)
x1 = np.linspace(np.min(h[:,1]-1), np.max(h[:,1]+1), 50)

X0, X1 = np.meshgrid(x0, x1)
Z = f(X0, X1)

plot3D(X0,X1,Z, initx, r.x, h)