# Numerical optimization

You will learn to solve non-convex multi-dimensional optimization problems using numerical optimization with multistart and nesting (**scipy.optimize**). You will learn simple function approximation using linear interpolation (**scipy.interp**). 

**Links:**

1. **scipy.optimize:** [overview](https://docs.scipy.org/doc/scipy/reference/optimize.html) + [tutorial](https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html)
2. **scipy.interp:** [overview](https://docs.scipy.org/doc/scipy/reference/interpolate.html) + [tutorial](https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html)

**Useful note:** [Numerical Optimization in MATLAB](http://www.google.com/url?q=http%3A%2F%2Fweb.econ.ku.dk%2Fmunk-nielsen%2Fnotes%2FnoteOptimization.pdf&sa=D&sntz=1&usg=AFQjCNHX4tHx2_YsNaIt5FB5MBU5cfcS8g) (by Anders Munk-Nielsen)


```python
from types import SimpleNamespace
import numpy as np
import scipy as sp
from scipy import linalg
from scipy import optimize
from scipy import interpolate
import sympy as sm

%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
```

# Introduction

All **optimization problems** are characterized by:

1. Control vector (choices), $\boldsymbol{x} \in \mathbb{R}^k$
2. Objective function (payoff) to minimize, $f:\mathbb{R}^k \rightarrow \mathbb{R}$ (differentiable or not)
3. Constraints, i.e. $\boldsymbol{x}  \in C \subseteq \mathbb{R}^k$ (linear or non-linear interdependence)

Note that $f$ might also take other inputs (parameters or a dataset), but these are fixed, and therefore not variables we optimize over.

**Maximization** is just **minimization** of $-f$. 

All **optimizers** (minimizers) follow the structure:

1. Make initial guess
2. Evaluate the function (and perhaps gradients)
3. Check for convergence
4. Update guess and return to step 2

**Convergence:** "Small" change in function value since last iteration (or "zero" gradient).

**Characteristics** of optimizers:

1. Use gradients or not.
2. Allow for specifying bounds.
3. Allow for specifying general constraints.

**Gradients** provide useful information, but can be costly to compute (using analytical formula or numerically).

**Penalty terms** can (sometimes) instead be used to enforce bounds and constraints.

**Optimizers** you should know:

1. **Nelder-Mead:** 
 * **Pro:** Robust (to e.g. noise in objective function) and does not require derivatives.
 * **Con:** Slow convergence. No bounds or constraints.
2. **Newton-CG:**
 * **Pro:** Require few iterations. Very precise with analytical hessian for smooth functions.
 * **Con:** Costly computation of hessian. No bounds or constraints.
3. **BFGS:** (like newton, but with smart computation of hessian)
  * **Pro:** Require few function evaluations. 
  * **Con:** No bounds or constraints.
4. **L-BFGS-B:** Like BFGS, but allows for bounds.
5. **SLSQP:**
  * **Pro:** Bounds and constraints in multiple dimensions.
  * **Con:** Not as efficient as BFGS.

## Gradient based optimizers

Let us look at the idea behind gradient based optimizers.

**One dimensional intuition:** Consider the second-order Taylor approximation around $x_n$:

$$ 
f_T(x) = f_T(x_n + \Delta x) \approx f(x_n)+ f^{\prime}(x_n) \Delta x + \frac{1}{2} f^{\prime\prime}(x_n) (\Delta x)^2
$$

Find the minimum wrt. to $\Delta x$ by solving the FOC:

$$
0 = \frac{d}{d\Delta x} f_T(x) = f^{\prime}(x_n) + f^{\prime\prime}(x_n) \Delta x \Leftrightarrow \Delta x = -\frac{f^{\prime}(x_n)}{f^{\prime\prime}(x_n)}
$$

**Algorithm:** `minimize_newton()`

1. Choose tolerance $\epsilon>0$, guess on $\boldsymbol{x}_0$, compute $f(\boldsymbol{x}_0)$, and set $n=1$.
2. Compute $\nabla f(\boldsymbol{x}_{n-1})$ (gradient/jacobian) and $\boldsymbol{H}f(\boldsymbol{x}_{n-1})$ (hessian).
3. Compute new guess

  $$ 
  \boldsymbol{x}_{n} = \boldsymbol{x}_{n-1} - [\boldsymbol{H}f(\boldsymbol{x}_{n-1})]^{-1} \nabla f(\boldsymbol{x}_{n-1})
  $$

3. If $|f(\boldsymbol{x}_n)-f(\boldsymbol{x}_{n-1})| < \epsilon$ then stop.
5. Set $n = n + 1$ and return to step 2.


```python
def minimize_newton(f,x0,jac,hess,max_iter=500,tol=1e-8):
    """ minimize function with Newtons' algorithm
        
    Args:

        f (callable): function
        x0 (np.ndarray): initial values
        jac (callable): jacobian
        hess (callable): hessian
        max_iter (int): maximum number of iterations
        tol (float): tolerance
        
    Returns:
    
        x (np.ndarray): minimum
        n (int): number of iterations used
        
    """
    
    # step 1: initialize
    x = x0
    fx = f(x0)
    n = 1
    
    # step 2-5: iteration
    while n < max_iter:
        
        x_prev = x
        fx_prev = fx
        
        # step 2: evaluate gradient and hessian
        jacx = jac(x_prev)
        hessx = hess(x_prev)
        
        # step 3: update x
        inv_hessx = linalg.inv(hessx)        
        x = x_prev - inv_hessx@jacx
     
        # step 4: check convergence
        fx = f(x)
        if abs(fx-fx_prev) < tol:
            break
            
        # step 5: increment n
        n += 1
        
    return x,n
```

**Algorithm:** `minimize_gradient_descent()`

1. Choose tolerance $\epsilon>0$, potential step sizes, $ \boldsymbol{\alpha} = [\alpha_0,\alpha_1,\dots,\alpha_\#]$, guess on $\boldsymbol{x}_0$, compute $f(\boldsymbol{x}_0)$, and set $n=1$.
2. Compute $\nabla f(\boldsymbol{x}_{n-1})$.
3. Find good step size:

  $$ 
  \alpha^{\ast} = \arg \min_{\alpha \in \boldsymbol{\alpha}}  f(\boldsymbol{x}_{n-1} - \alpha \nabla f(\boldsymbol{x}_{n-1}))
  $$

4. Compute new guess:

  $$
  \boldsymbol{x}_{n} = \boldsymbol{x}_{n-1} - \alpha^{\ast} \nabla f(\boldsymbol{x}_{n-1})
  $$

5. If $|f(\boldsymbol{x}_n)-f(\boldsymbol{x}_{n-1})| < \epsilon$ then stop.
6. Set $n = n + 1$ and return to step 2.


```python
def minimize_gradient_descent(f,x0,jac,alphas=[0.01,0.05,0.1,0.25,0.5,1],max_iter=500,tol=1e-8):
    """ minimize function with gradient descent
        
    Args:

        f (callable): function
        x0 (np.ndarray): initial values
        jac (callable): jacobian
        alpha (list): potential step sizes
        max_iter (int): maximum number of iterations
        tol (float): tolerance
        
    Returns:
    
        x (np.ndarray): minimum
        n (int): number of iterations used
        
    """
    
    # step 1: initialize
    x = x0
    fx = f(x0)
    n = 1
    
    # step 2-6: iteration
    while n < max_iter:
            
        x_prev = x
        fx_prev = fx
        
        # step 2: evaluate gradient
        jacx = jac(x)
        
        # step 3: find good step size (line search)
        fx_ast = np.inf
        alpha_ast = np.nan
        for alpha in alphas:
            x = x_prev - alpha*jacx
            fx = f(x)
            if fx < fx_ast:
                fx_ast = fx
                alpha_ast = alpha
        
        # step 4: update guess
        x = x_prev - alpha_ast*jacx
                            
        # step 5: check convergence
        fx = f(x)
        if abs(fx-fx_prev) < tol:
            break
            
        # d. update i
        n += 1
        
    return x,n
```

**Many generalizations:**

1. Use both Hessian and line search
2. Stop line search when improvement found
3. Limit attention to a "trust-region"

etc. etc. etc. etc.

# Example: The rosenbrock function

Consider the **rosenbrock function**:

$$ 
f(\boldsymbol{x}) = f(x_1,x_2) =0.5(1-x_{1})^{2}+(x_{2}-x_{1}^{2})^{2}
$$

with **jacobian** (gradient)

$$ 
\nabla f(\boldsymbol{x})=\begin{bmatrix}\frac{\partial f}{\partial x_{1}}\\
\frac{\partial f}{\partial x_{2}}
\end{bmatrix}=\begin{bmatrix}-(1-x_{1})-4x_{1}(x_{2}-x_{1}^{2})\\
2(x_{2}-x_{1}^{2})
\end{bmatrix}
$$

and **hessian**:

$$
\boldsymbol{H}f(\boldsymbol{x})=\begin{bmatrix}\frac{\partial f}{\partial x_{1}x_{1}} & \frac{\partial f}{\partial x_{1}x_{2}}\\
\frac{\partial f}{\partial x_{1}x_{2}} & \frac{\partial f}{\partial x_{2}x_{2}}
\end{bmatrix}=\begin{bmatrix}1-4x_{2}+12x_{1}^{2} & -4x_{1}\\
-4x_{1} & 2
\end{bmatrix}
$$

**Note:** Minimum is at $(1,1)$ where $f(1,1)=0$.

**Check jacobian and hessian:**


```python
x1 = sm.symbols('x_1')
x2 = sm.symbols('x_2')
f = 0.5*(1.0-x1)**2 + (x2-x1**2)**2
```


```python
Df = sm.Matrix([sm.diff(f,i) for i in [x1,x2]])
Df
```




$\displaystyle \left[\begin{matrix}- 4 x_{1} \left(- x_{1}^{2} + x_{2}\right) + 1.0 x_{1} - 1.0\\- 2 x_{1}^{2} + 2 x_{2}\end{matrix}\right]$




```python
Hf = sm.Matrix([[sm.diff(f,i,j) for j in [x1,x2]] for i in [x1,x2]])
Hf
```




$\displaystyle \left[\begin{matrix}12 x_{1}^{2} - 4 x_{2} + 1.0 & - 4 x_{1}\\- 4 x_{1} & 2\end{matrix}\right]$



**Implementation:**


```python
def _rosen(x1,x2):
    return 0.5*(1.0-x1)**2+(x2-x1**2)**2
def rosen(x):
    return _rosen(x[0],x[1])
def rosen_jac(x):
    return np.array([-(1.0-x[0])-4*x[0]*(x[1]-x[0]**2),2*(x[1]-x[0]**2)])
def rosen_hess(x):
    return np.array([[1-4*x[1]+12*x[0]**2,-4*x[0]],[-4*x[0],2]])
```

**3D Plot:**


```python
# a. grids
x1_vec = np.linspace(-2,2,500)
x2_vec = np.linspace(-2,2,500)
x1_grid,x2_grid = np.meshgrid(x1_vec,x2_vec,indexing='ij')
rosen_grid = _rosen(x1_grid,x2_grid)

# b. main
fig = plt.figure()
ax = fig.add_subplot(1,1,1,projection='3d')
cs = ax.plot_surface(x1_grid,x2_grid,rosen_grid,cmap=cm.jet)

# c. add labels
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$u$')

# d. invert xaxis
ax.invert_xaxis()

# e. add colorbar
fig.colorbar(cs);
```

**Contour plot:**


```python
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
levels = [1e-6,5*1e-6,1e-5,5*1e-5,1e-4,5*1e-4,1e-3,5*1e-3,1e-2,5*1e-2,1,2,4,6,8,12,16,20]
cs = ax.contour(x1_grid,x2_grid,rosen_grid,levels=levels,cmap=cm.jet)
fig.colorbar(cs);
```

**Newton:**


```python
x0 = np.array([5,4])
x,n = minimize_newton(rosen,x0,rosen_jac,rosen_hess)
print(n,x,rosen(x))
```

    6 [1. 1.] 0.0


**Gradient descent:**


```python
x0 = np.array([5,4])
x,n = minimize_gradient_descent(rosen,x0,rosen_jac,alphas=[0.01,0.05,0.1,0.25,0.5,1])
print(n,x,rosen(x))
```

    173 [1.00020519 1.00053964] 3.7750814497569406e-08


**Questions:** Any ideas for getting the gradient descent optimizer to converge faster?

## Scipy minimizers

**Preperation I:** Function for collecting infomation while running optimizing:


```python
# complicated -> not necessary to understand it
def collect(x):
    
    # globals used to keep track across iterations
    global evals # set evals = 0 before calling optimizer
    global x0
    global x1s
    global x2s
    global fs
    
    # a. initialize list
    if evals == 0:
        x1s = [x0[0]] 
        x2s = [x0[1]]
        fs = [rosen(x0)]
        
    # b. append trial values
    x1s.append(x[0])
    x2s.append(x[1])
    fs.append(rosen(x))
    
    # c. increment number of evaluations
    evals += 1
    
```

**Preperation II:** Function plotting the collected information:


```python
# complicated -> not necessary to understand it
def contour():
    
    global evals
    global x1s
    global x2s
    global fs
    
    # a. contour plot
    fig = plt.figure(figsize=(10,4))
    ax = fig.add_subplot(1,2,1)
    levels = [1e-6,5*1e-6,1e-5,5*1e-5,1e-4,5*1e-4,1e-3,5*1e-3,1e-2,5*1e-2,1,2,4,6,8,12,16,20]
    cs = ax.contour(x1_grid,x2_grid,rosen_grid,levels=levels,cmap=cm.jet)
    fig.colorbar(cs)
    ax.plot(x1s,x2s,'-o',ms=4,color='black')
    ax.set_xlabel('$x_1$')
    ax.set_ylabel('$x_2$')
    
    # b. function value
    ax = fig.add_subplot(1,2,2)
    ax.plot(np.arange(evals+1),fs,'-o',ms=4,color='black')
    ax.set_xlabel('iteration')
    ax.set_ylabel('function value')
    
```

**Nelder-Mead**


```python
evals = 0 # global used in "collect"
x0 = [-1.5,-1]
result = optimize.minimize(rosen,x0,
                           method='Nelder-Mead',
                           callback=collect, # call collect() before each iteration
                           options={'disp':True}) # display the results
contour()
```

> **Note:** Does not require a gradient. Slow convergence close to target.
>
> **Iterations:** How many steps the algorithm has taken.
>
> **Function evaluations:** Will be higher than iterations. Used to compute next step.

We can also **print the information on results:**


```python
print(result)
```

     final_simplex: (array([[0.9999932 , 0.99998507],
           [1.00000809, 1.00003102],
           [1.00003176, 1.00007406]]), array([2.48761443e-11, 2.53004854e-10, 6.15284951e-10]))
               fun: 2.4876144250933352e-11
           message: 'Optimization terminated successfully.'
              nfev: 105
               nit: 57
            status: 0
           success: True
                 x: array([0.9999932 , 0.99998507])


We can also acess specific information of the result object: 


```python
result.nit
```




    57



**Newton** (with analytical hessian)


```python
evals = 0 # global used in "collect"
x0 = [-1.5,-1]
result = optimize.minimize(rosen,x0,jac=rosen_jac,hess=rosen_hess,
                           method='Newton-CG',
                           callback=collect,
                           options={'disp':True})
contour()
```

> **Note:** Smoother and faster.

**Newton** (with numerical hessian computed by scipy)


```python
evals = 0 # global used in "collect"
x0 = [-1.5,-1]
result = optimize.minimize(rosen,x0,jac=rosen_jac,
                           method='Newton-CG',
                           callback=collect,
                           options={'disp':True})
contour()
```

> **Note:** Same as above, but gradient evaluations instead of hessian evaluations.

**BFGS** (with analytical gradient)


```python
evals = 0 # global used in "collect"
x0 = [-1.5,-1]
result = optimize.minimize(rosen,x0,jac=rosen_jac,
                           method='BFGS',
                           callback=collect,
                           options={'disp':True})
contour()
```

> **Note:** Non-smooth, but fast. Very low number of function evaluations.

**BFGS** (with numerical gradient computed by scipy)


```python
evals = 0 # global used in "collect"
x0 = [-1.5,-1]
result = optimize.minimize(rosen,x0, # no jac= specified
                           method='BFGS',
                           callback=collect,
                           options={'disp':True})
contour()
```

> **Note:** Same as above, but more function evaluations.

**L-BFGS-B** (with analytical gradient)


```python
evals = 0 # global used in "collect"
x0 = [-1.5,-1]
result = optimize.minimize(rosen,x0,jac=rosen_jac,
                           method='L-BFGS-B',
                           bounds=((-3,3),(-3,3)),
                           callback=collect,
                           options={'disp':True})
contour()
```

**SLSQP**


```python
evals = 0 # global used in "collect"
x0 = [-1.5,-1]
result = optimize.minimize(rosen,x0,jac=rosen_jac,
                           method='SLSQP',
                           bounds=((-2,2),(-2,2)),
                           callback=collect,
                           options={'disp':True})
contour()
```

## Controling the optimizers

> **Note:** See the settings for each optimizer in the [documention](https://docs.scipy.org/doc/scipy/reference/optimize.html).

We can lower the **tolerance**:


```python
evals = 0
x0 = [-1.5,-1]
result = optimize.minimize(rosen,x0,
                           method='BFGS',
                           callback=collect,
                           options={'disp':True,'gtol':1e-8}) # note this
contour()
```

We can change the **maximum number of iterations**:


```python
evals = 0
x0 = [-1.5,-1]
result = optimize.minimize(rosen,x0,
                           method='BFGS',
                           callback=collect,
                           options={'disp':True,'maxiter':5}) # note this and warning
contour()
```

**Question:** Can we make the program stop if the maximum number of iterations is too low?

# Sombrero function: Local minima and multistart

Consider the **sombrero** function

$$
f(x_1,x_2) = g\Big(\sqrt{x_1^2 + x_2^2}\Big)
$$

where

$$
g(r) = -\frac{\sin(r)}{r+10^{-4}} + 10^{-4}r^2
$$

The **global minimum** of this function is (0,0). But the function also have (infinitely many) **local minima**. How to avoid these?


```python
def _sombrero(x1,x2):
    r = np.sqrt(x1**2 + x2**2)
    return -np.sin(r)/(r+1e-4) + 1e-4*r**2
    
sombrero = lambda x: _sombrero(x[0],x[1])
```

## 3D plot


```python
# a. grids
x1_vec = np.linspace(-15,15,500)
x2_vec = np.linspace(-15,15,500)
x1_grid_sombrero,x2_grid_sombrero = np.meshgrid(x1_vec,x2_vec,indexing='ij')
sombrero_grid = _sombrero(x1_grid_sombrero,x2_grid_sombrero)

# b. main
fig = plt.figure()
ax = fig.add_subplot(1,1,1,projection='3d')
cs = ax.plot_surface(x1_grid_sombrero,x2_grid_sombrero,sombrero_grid,cmap=cm.jet)

# c. add labels
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$u$')

# d. invert xaxis
ax.invert_xaxis()

# e. colorbar
fig.colorbar(cs);
```

## Multi-start - BFGS

**Multi-start:** Draw many random starting values:


```python
np.random.seed(1986)
x0s = -15 + 30*np.random.uniform(size=(5000,2)) # in [-15,15]
xs = np.empty((5000,2))
fs = np.empty(5000)
```

Try to solve with **BFGS** starting from each of these:


```python
fopt = np.inf
xopt = np.nan
for i,x0 in enumerate(x0s):
    
    # a. optimize
    result = optimize.minimize(sombrero,x0,method='BFGS')
    xs[i,:] = result.x
    f = result.fun
    
    # b. print first 10 or if better than seen yet
    if i < 10 or f < fopt: # plot 10 first or if improving
        if f < fopt:
            fopt = f
            xopt = xs[i,:]
            
        print(f'{i:4d}: x0 = ({x0[0]:6.2f},{x0[1]:6.2f})',end='')
        print(f' -> converged at ({xs[i][0]:6.2f},{xs[i][1]:6.2f}) with f = {f:.14f}')
        
# best solution
print(f'\nbest solution:\n x = ({xopt[0]:6.2f},{xopt[1]:6.2f}) -> f = {fopt:.14f}')
```

       0: x0 = (  2.07,-12.68) -> converged at (  2.26,-13.84) with f = -0.05118272212318
       1: x0 = (  3.25, 11.92) -> converged at (  3.69, 13.53) with f = -0.05118272211067
       2: x0 = (  1.35,  6.08) -> converged at (  1.67,  7.53) with f = -0.12241421149687
       3: x0 = ( -3.42,  4.55) -> converged at ( -4.63,  6.17) with f = -0.12241421149362
       4: x0 = (  5.70,  6.56) -> converged at (  5.06,  5.82) with f = -0.12241421149696
       5: x0 = (  5.71,  9.20) -> converged at (  4.06,  6.56) with f = -0.12241421146807
       6: x0 = ( -6.49,  7.73) -> converged at ( -4.96,  5.91) with f = -0.12241421149702
       7: x0 = ( -5.77, -1.86) -> converged at ( -7.34, -2.37) with f = -0.12241421135584
       8: x0 = (  1.58, -2.01) -> converged at (  0.04, -0.05) with f = -0.99776236017000
       9: x0 = (  2.33, 11.77) -> converged at (  2.72, 13.76) with f = -0.05118272212317
      25: x0 = (  1.41, -3.37) -> converged at (  0.03, -0.06) with f = -0.99776236017097
     150: x0 = (-14.53,  8.53) -> converged at ( -0.06,  0.03) with f = -0.99776236017097
     223: x0 = (  2.87, -0.21) -> converged at ( -0.07,  0.00) with f = -0.99776236017098
     279: x0 = ( -6.90, -8.37) -> converged at ( -0.04, -0.05) with f = -0.99776236017098
     354: x0 = ( -1.24, -2.85) -> converged at ( -0.03, -0.06) with f = -0.99776236017098
    1588: x0 = (-12.93,-11.42) -> converged at ( -0.05, -0.04) with f = -0.99776236017098
    4700: x0 = (  1.06,  0.65) -> converged at (  0.06,  0.04) with f = -0.99776236017098
    
    best solution:
     x = (  0.06,  0.04) -> f = -0.99776236017098


The solver, wrongly, **converges to many of the local minima**:


```python
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.scatter(xs[:,0],xs[:,1]);
```

## Multi-start - Nelder-Mead

Try to solve with **Nelder-Mead** starting from each of these:


```python
fopt = np.inf
xopt = np.nan
for i,x0 in enumerate(x0s):
    
    # a. optimize
    result = optimize.minimize(sombrero,x0,method='Nelder-Mead')
    xs[i,:] = result.x
    f = result.fun
    
    # b. print first 10 or if better than seen yet
    if i < 10 or f < fopt: # plot 10 first or if improving
        if f < fopt:
            fopt = f
            xopt = xs[i,:]
        print(f'{i:4d}: x0 = ({x0[0]:6.2f},{x0[1]:6.2f})',end='')
        print(f' -> converged at ({xs[i][0]:6.2f},{xs[i][1]:6.2f}) with f = {f:.12f}')

# best solution
print(f'\nbest solution:\n x = ({xopt[0]:6.2f},{xopt[1]:6.2f}) -> f = {fopt:.12f}')
```

       0: x0 = (  2.07,-12.68) -> converged at (  2.20,-13.85) with f = -0.051182722123
       1: x0 = (  3.25, 11.92) -> converged at (  3.43, 13.60) with f = -0.051182722123
       2: x0 = (  1.35,  6.08) -> converged at (  1.39,  7.59) with f = -0.122414211497
       3: x0 = ( -3.42,  4.55) -> converged at ( -3.71,  6.76) with f = -0.122414211497
       4: x0 = (  5.70,  6.56) -> converged at (  6.69,  3.84) with f = -0.122414211497
       5: x0 = (  5.71,  9.20) -> converged at (  4.98, 13.11) with f = -0.051182722123
       6: x0 = ( -6.49,  7.73) -> converged at ( -7.34,  2.37) with f = -0.122414211497
       7: x0 = ( -5.77, -1.86) -> converged at ( -7.45, -2.01) with f = -0.122414211497
       8: x0 = (  1.58, -2.01) -> converged at ( -0.02,  0.06) with f = -0.997762360171
       9: x0 = (  2.33, 11.77) -> converged at (  2.36, 13.83) with f = -0.051182722123
    
    best solution:
     x = ( -0.02,  0.06) -> f = -0.997762360171



```python
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.scatter(xs[:,0],xs[:,1]);
```

## Is there a better solution than multi-start?

**In short:** No.

**Potential improvement:** Use information from previous run to determine, where to look next. Fundamental trade-off between:

1. **Exploitation.** Focus on areas where previous evaluations returned low function values.
2. **Exploration.** Focus on completely new areas. 

**Heuristic:** If the same optimum is obtained for many starting values, this is a good sign for it being the global optimum.

**Further discussion**: [Benchmarking Global Optimizers](https://fguvenendotcom.files.wordpress.com/2019/09/agk2019-september-nber-submit.pdf) ([code](https://github.com/serdarozkan/TikTak#tiktak))

# Constraints

## In general

Consider the **constrained problem**:

$$
\min_{x_1,x_2,x_3,x_4} x_1x_4(x_1+x_2+x_3) + x_3
$$

subject to

$$
\begin{aligned}
x_1x_2x_3x_4 &\geq 25 \\
x_1^2+x_2^2+x_3^2+x_4^2 &= 40 \\
1 \leq x_1,x_2,x_3,x_4 &\leq 5
\end{aligned}
$$

Define **objective** and **constraints**:


```python
def _objective(x1,x2,x3,x4):
    return x1*x4*(x1+x2+x3)+x3

def objective(x):
    return _objective(x[0],x[1],x[2],x[3])

def ineq_constraint(x):
    return x[0]*x[1]*x[2]*x[3]-25.0 # violated if negative

def eq_constraint(x):
    sum_eq = 40.0
    for i in range(4):
        sum_eq = sum_eq - x[i]**2
    return sum_eq # must equal zero
```


```python
# a. setup
bound = (1.0,5.0)
bounds = (bound, bound, bound, bound)
ineq_con = {'type': 'ineq', 'fun': ineq_constraint} 
eq_con = {'type': 'eq', 'fun': eq_constraint}

# b. call optimizer
x0 = (40**(1/8),40**(1/8),40**(1/8),40**(1/8)) # fit the equality constraint
result = optimize.minimize(objective,x0,
                             method='SLSQP',
                             bounds=bounds,
                             constraints=[ineq_con,eq_con],
                             options={'disp':True})

print('\nx = ',result.x)
```

    Optimization terminated successfully    (Exit mode 0)
                Current function value: 17.0140172890441
                Iterations: 9
                Function evaluations: 46
                Gradient evaluations: 9
    
    x =  [1.         4.74299964 3.82114999 1.37940829]


**Alternative:** Extend the **objective function with a penalty term**, where guesses outside the allowed bounds and constraints are projected into the allowed region, but a (large) penalty term is added to discourage this. Solve this problem with an unconstrained solver.

## Economic application

Consider the following **consumption-saving problem**:

$$
\begin{aligned}
V(a_0) &= \max_{c_1,c_2,c_3} \frac{c_{1}^{1-\rho}}{1-\rho} + \beta \frac{c_{2}^{1-\rho}}{1-\rho} + \beta^2\frac{c_{3}^{1-\rho}}{1-\rho} + \beta^2\nu\frac{(a_{3}+\kappa)^{1-\rho}}{1-\rho}\\
\text{s.t.} \\
&\text{s.t.}&\\
m_1 &= (1+r)a_0 + y_1\\
a_1 &= m_1-c_1\\
m_2 &= (1+r)a_1 + y_2\\
a_2 &= m_2-c_2\\
m_3 &= (1+r)a_2 + y_3\\
a_3 &= m_3-c_3\\
c_1,c_2,c_3 &\geq 0\\
a_1,a_2,a_3 &\geq 0\\
\end{aligned}
$$

where 

* $m_t$ is cash-on-hand in period $t\in\{1,2,\dots,T\}$
* $c_t$ is consumption $t$
* $a_t$ is end-of-period assets and income in period $t$
* ${y_t}$ is income in period $t$
* $\beta > 0$ is the discount factor
* $r > -1$ is the interest rate 
* $\rho > 1$ is the CRRA coefficient
* $\nu > 0 $ is the strength of the bequest motive
* $\kappa > 0$ is the degree of luxuriousness in the bequest motive  
* $a_t\geq0$ is a no-borrowing constraint.


**Guide to solve such problem:**

1. Setup parameters
2. Formulate objective function
3. Determine how to handle constraints
4. Call optimizer

**Parameters:**


```python
par = SimpleNamespace()
par.a0 = 0.5
par.beta = 0.94
par.r = 0.04
par.rho = 8
par.kappa = 0.5
par.nu = 0.1
par.y = [1,2,4]
par.T = 3
```

**Objetive function:**


```python
def obj(c,par,full_return=False):
    
    # objective function with penalty term
    
    # a. allocate
    a = np.zeros(par.T) # end-of-period assets
    m = np.zeros(par.T) # cash-on-hand
    cb = np.zeros(par.T) # bounded
    
    # b. bound consumption and penalty
    penalty = 0.0
    for t in range(par.T):
        
        # i. lagged assets
        a_lag = a[t-1] if t > 0 else par.a0
        
        # ii. cash-on-hand
        m[t] = (1+par.r)*a_lag + par.y[t]
        
        # ii. bounded consumption
        if c[t] < 0:
            penalty += 10_000*np.abs(c[t]-0.0)
            cb[t] = 0
        elif c[t] > m[t]:
            penalty += 10_000*np.abs(c[t]-m[t])
            cb[t] = m[t]
        else:
            cb[t] = c[t]
        
        # d. end-of-period assets 
        a[t] = m[t] - cb[t]
            
    # c. utility
    utility = 0.0
    
    # i. consumption
    for t in range(par.T):
        utility += par.beta**t*(cb[t]**(1-par.rho))/(1-par.rho)
    
    # ii. bequest
    utility += par.beta**(par.T-1)*par.nu*(a[-1]+par.kappa)**(1-par.rho)/(1-par.rho)
        
    # d. return negative utility + penalty
    if full_return:
        return utility,m,a
    else:
        return -utility + penalty
```

**Solve:**


```python
def solve(par):
    
    # a. initial geuss
    x0 = [par.a0/par.T,par.a0/par.T,par.a0/par.T]
    
    # b. solve
    results = optimize.minimize(obj,x0,args=(par,),method='nelder-mead')
    assert results.success
    print(f'solved in {results.nit} iteratoons [{results.nfev} function evaluations]')
    
    # c. details
    c = results.x
    utility,m,a = obj(c,par,full_return=True)
    print(f't = 0: a = {par.a0:.4f}')
    for t in range(par.T):
        print(f't = {t+1}: y = {par.y[t]:.4f}, m = {m[t]:.4f}, c = {c[t]:.4f}, a = {a[t]:.4f}')    
    print(f'utility = {utility:.8f}')
```


```python
solve(par)
```

    solved in 157 iteratoons [287 function evaluations]
    t = 0: a = 0.5000
    t = 1: y = 1.0000, m = 1.5200, c = 1.5200, a = 0.0000
    t = 2: y = 2.0000, m = 2.0000, c = 2.0000, a = 0.0000
    t = 3: y = 4.0000, m = 4.0000, c = 2.5675, a = 1.4325
    utility = -0.00896688


**What happens if the income path is reversed?**


```python
par.y = list(reversed(par.y))
solve(par)
```

    solved in 138 iteratoons [242 function evaluations]
    t = 0: a = 0.5000
    t = 1: y = 4.0000, m = 4.5200, c = 2.1948, a = 2.3252
    t = 2: y = 2.0000, m = 4.4182, c = 2.1885, a = 2.2297
    t = 3: y = 1.0000, m = 3.3189, c = 2.1823, a = 1.1365
    utility = -0.00207765


**Question I:** Could we easily allowing for borrowing, i.e.

$$
\begin{aligned}
V(a_0) &= \max_{c_1,c_2,c_3} \frac{c_{1}^{1-\rho}}{1-\rho} + \beta \frac{c_{2}^{1-\rho}}{1-\rho} + \beta^2\frac{c_{3}^{1-\rho}}{1-\rho} + \beta^2\nu\frac{(a_{3}+\kappa)^{1-\rho}}{1-\rho}\\
\text{s.t.} \\
&\text{s.t.}&\\
m_1 &= (1+r)a_0 + y_1\\
a_1 &= m_1-c_1\\
m_2 &= (1+r)a_1 + y_2\\
a_2 &= m_2-c_2\\
m_3 &= (1+r)a_2 + y_3\\
a_3 &= m_3-c_3\\
c_1,c_2,c_3 &\geq 0\\
a_1,a_2 &\geq -\lambda\\
a_3 &\geq 0
\end{aligned}
$$

**Question II:** Could we easily extend the problem to more periods?

$$
\begin{aligned}
V(a_0) &= \max_{c_1,c_2,\dots c_T} \sum_{t=1}^T \beta^{t-1} \frac{c_{t}^{1-\rho}}{1-\rho} + \beta^{T+1}\nu\frac{(a_{T}+\kappa)^{1-\rho}}{1-\rho}\\
\text{s.t.} \\
&\text{s.t.}&\\
m_t &= (1+r)a_{t-1} + y_t\\
c_t &\geq 0\\
a_t &\geq 0
\end{aligned}
$$

**Follow-up question:** What is the problem for $T \rightarrow \infty$?

# Interpolation

**Intermezzo:** To consider dynamic optimization problems, we need to think about interpolation.

**Inputs:**

1. Sorted vector of known points (grid vector), $G = \{G_i\}_{i=0}^{n-1}$
2. Vector of known values (at these points), $F = \{F_i = f(G_i)\}_{i=0}^{n-1}$
3. A new point, `x`

**Algorithm:** `linear_interpolate()`
1. Determine `i`  such that

$$
G_i \leq x < G_{i+1}
$$

2. Compute interpolated value by

$$
y =  F_{i} + \frac{F_{i+1}-F_{i}}{G_{i+1}-G_{i}}(x-G_{i})
$$

**Extrapolation:**

1. Below where $x < G_1$: 

$$
y =  F_{0} + \frac{F_{1}-F_{0}}{G_{1}-G_{0}}(x-G_{0})
$$

2. Above where $x > G_{n-2}$: 

$$
y =  F_{n-2} + \frac{F_{n-1}-F_{n-2}}{G_{n-1}-G_{n-2}}(x-G_{n-2})
$$


```python
def linear_interpolate(G,F,x):
    """ linear interpolation (and extrapolation)
    
    Args:
    
        G (np.ndarray): known points
        F (np.ndarray): known values
        x (float): point to be interpolated
        
    Returns:
    
        y (float): intepolated value
    
    """
    
    assert len(G) == len(F)
    n = len(G)
    
    # a. find index in known points
    if x < G[1]: # exprapolation below
        i = 0
    elif x > G[-2]: # extrapolation above
        i = n-2
    else: # true interpolation
        
        # search
        i = 0 
        while x >= G[i+1] and i < n-1:
            i += 1
        
        assert x >= G[i]
        assert x < G[i+1]

    # b. interpolate
    diff_G = G[i+1]-G[i]
    diff_F = F[i+1]-F[i]
    slope = diff_F/diff_G
    y = F[i] + slope*(x-G[i])
    
    return y
```

## Example

Consider the following function and known points:


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

G = np.linspace(-5,10,6)
F = f(G)
```

**Simple test:**


```python
for x in [-2.3,4.1,7.5,9.1]:
    true = f(x)
    y = linear_interpolate(G,F,x)
    print(f'x = {x:4.1f} -> true = {true:6.1f}, interpolated = {y:6.1f}')
```

    x = -2.3 -> true = -176.2, interpolated = -193.5
    x =  4.1 -> true =  -28.6, interpolated =  -27.7
    x =  7.5 -> true =  -40.1, interpolated =  -24.5
    x =  9.1 -> true =   24.1, interpolated =   50.7


**Scipy.interpolate:** Use the *RegularGridInterpolator*  


```python
# a. construct interpolation function
interp_func = interpolate.RegularGridInterpolator([G],F,
                                                  bounds_error=False,
                                                  fill_value=None)

# bounds_error=False and fill_value=None allow for extrapolation

# b. interpolate
grid = np.linspace(-7,12,500)
interp_values = interp_func(grid)

# c. evaluate true values
true_values = f(grid)

# d. plot true and interpolated values
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(G,F,'o',label='known points')
ax.plot(grid,true_values,'-',lw=1,label='true function')
ax.plot(grid,interp_values,'-',lw=1,label='interpolated values')
ax.legend(loc='lower right',facecolor='white',frameon=True);
```

**Note:**

1. Linear interpolation works best when the function does not curve too much.
2. Extrapolation is much worse than interpolation.

**Multiple dimensions:** Same principle, ``interpolate.RegularGridInterpolator([G1,G2,G3],F)``.

# Dynamic optimization problems

The following subject is hard. But also extremely useful. *If you master this, you can solve (almost) all economic models you meet on your way in life*.

## Problem formulation

Consider a **household** living in two periods.

In the **second period** it gets utility from **consuming** and **leaving a bequest** (warm glow),

$$
\begin{aligned}
v_{2}(m_{2})&= \max_{c_{2}}\frac{c_{2}^{1-\rho}}{1-\rho}+\nu\frac{(a_2+\kappa)^{1-\rho}}{1-\rho}\\
\text{s.t.} \\
a_2 &= m_2-c_2 \\
a_2 &\geq 0
\end{aligned}
$$

where 

* $m_2$ is cash-on-hand 
* $c_2$ is consumption
* $a_2$ is end-of-period assets 
* $\rho > 1$ is the risk aversion coefficient
* $\nu > 0 $ is the strength of the bequest motive
* $\kappa > 0$ is the degree of luxuriousness in the bequest motive  
* $a_2\geq0$ ensures the household *cannot* die in debt

The **value function** $v(m_2)$ measures the household's value of having $m_2$ at the beginning of period 2.


```python
def utility(c,par):
    return c**(1-par.rho)/(1-par.rho)

def bequest(m,c,par):
    return par.nu*(m-c+par.kappa)**(1-par.rho)/(1-par.rho)

def v2(c2,m2,par):
    return utility(c2,par) + bequest(m2,c2,par)
```

In the **first period**, the household gets utility from consuming and takes into account that it will also live in the next-period, where it receives a stochastic income,

$$
\begin{aligned}
v_1(m_1)&=\max_{c_1}\frac{c_{1}^{1-\rho}}{1-\rho}+\beta\mathbb{E}_{1}\left[v_2(m_2)\right]\\&\text{s.t.}&\\
a_1&=m_1-c_1\\
m_2&= (1+r)(m_1-c_1)+y_2 \\
y_{2}&= \begin{cases}
1-\Delta & \text{with prob. }0.5\\
1+\Delta & \text{with prob. }0.5 
\end{cases}\\
a_1&\geq0
\end{aligned}
$$

where

* $m_1$ is cash-on-hand in period 1
* $c_1$ is consumption in period 1
* $a_1$ is end-of-period assets in period 1
* $\beta > 0$ is the discount factor
* $\mathbb{E}_1$ is the expectation operator conditional on information in period 1
* $y_2$ is income in period 2
* $\Delta \in (0,1)$ is the level of income risk (mean-preserving)
* $r$ is the interest rate
* $a_1\geq0$ ensures the household *cannot* borrow


```python
def v1(c1,m1,par,v2_interp):
    
    # a. v2 value, if low income
    m2_low = (1+par.r)*(m1-c1) + 1-par.Delta
    v2_low = v2_interp([m2_low])[0]
    
    # b. v2 value, if high income
    m2_high = (1+par.r)*(m1-c1) + 1+par.Delta
    v2_high = v2_interp([m2_high])[0]
    
    # c. expected v2 value
    prob_low = 0.5
    prob_high = 0.05
    expected_v2 = prob_low*v2_low + prob_high*v2_high
    
    # d. total value
    return utility(c1,par) + par.beta*expected_v2
```

## Solve household problem

Choose **parameters**:


```python
par = SimpleNamespace()
par.rho = 8
par.kappa = 0.5
par.nu = 0.1
par.r = 0.04
par.beta = 0.94
par.Delta = 0.5
```

**Solve second period:**


```python
def solve_period_2(par):

    # a. grids
    m2_vec = np.linspace(1e-4,5,500)
    v2_vec = np.empty(500)
    c2_vec = np.empty(500)

    # b. solve for each m2 in grid
    for i,m2 in enumerate(m2_vec):

        # i. objective
        obj = lambda x: -v2(x[0],m2,par)

        # ii. initial value (consume half)
        x0 = m2/2

        # iii. optimizer
        result = optimize.minimize(obj,[x0],method='L-BFGS-B',bounds=((1e-8,m2),))

        # iv. save
        v2_vec[i] = -result.fun
        c2_vec[i] = result.x
        
    return m2_vec,v2_vec,c2_vec

# solve
m2_vec,v2_vec,c2_vec = solve_period_2(par)

# illustration
fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(1,2,1)
ax.plot(m2_vec,c2_vec)
ax.set_xlabel('$m_2$')
ax.set_ylabel('$c_2$')
ax.set_title('consumption function in period 2')

ax = fig.add_subplot(1,2,2)
ax.plot(m2_vec,v2_vec)
ax.set_xlabel('$m_2$')
ax.set_ylabel('$v_2$')
ax.set_title('value function in period 2')
ax.set_ylim([-40,1]);
```

**Note:** We now solve for the consumption function, rather than a specific optimum.

**Question:** Why is there a kink in the consumption function?

**Construct interpolator:**


```python
v2_interp = interpolate.RegularGridInterpolator([m2_vec], v2_vec,
                                                bounds_error=False,fill_value=None)
```

**Solve first period:**


```python
def solve_period_1(par,v2_interp):

    # a. grids
    m1_vec = np.linspace(1e-8,4,100)
    v1_vec = np.empty(100)
    c1_vec = np.empty(100)
    
    # b. solve for each m1 in grid
    for i,m1 in enumerate(m1_vec):
        
        # i. objective
        obj = lambda x: -v1(x[0],m1,par,v2_interp)
        
        # ii. initial guess (consume half)
        x0 = m1/2
        
        # iii. optimize
        result = optimize.minimize(obj,[x0],method='L-BFGS-B',bounds=((1e-12,m1),))
        
        # iv. save
        v1_vec[i] = -result.fun
        c1_vec[i] = result.x[0]
     
    return m1_vec,v1_vec,c1_vec

# solve
m1_vec,v1_vec,c1_vec = solve_period_1(par,v2_interp)

# illustrate
fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(1,2,1)
ax.plot(m1_vec,c1_vec)
ax.set_xlabel('$m_1$')
ax.set_ylabel('$c_1$')
ax.set_title('consumption function in period 1')

ax = fig.add_subplot(1,2,2)
ax.plot(m1_vec,v1_vec)
ax.set_xlabel('$m_1$')
ax.set_ylabel('$c_1$')
ax.set_title('value function in period 1')
ax.set_ylim([-40,1]);
```

**Summary:** We can summarize what we have done in a single function doing:

1. Solve period 2 (i.e. find $v_2(m_2)$ og $c_2(m_2)$)
2. Construct interpolator of $v_2(m_2)$
3. Solve period 1 (i.e. find $v_1(m_1)$ og $c_1(m_1)$)


```python
def solve(par):
    
    # a. solve period 2
    m2_vec,v2_vec,c2_vec = solve_period_2(par)
    
    # b. construct interpolator
    v2_interp = interpolate.RegularGridInterpolator([m2_vec], v2_vec,
        bounds_error=False,fill_value=None)
    
    # b. solve period 1
    m1_vec,v1_vec,c1_vec = solve_period_1(par,v2_interp)
    
    return m1_vec,c1_vec,m2_vec,c2_vec
```

**Plot consumption function for various level of income risk**, i.e varios $\Delta$


```python
fig = plt.figure()
ax = fig.add_subplot(1,1,1)

_Delta = par.Delta
for Delta in [0.05,0.15,0.25]:
    par.Delta = Delta
    m1_vec,c1_vec,m2_vec,c2_vec = solve(par)
    ax.plot(m1_vec,c1_vec,label=f'$\Delta = {Delta}$')

# reset
par.Delta = _Delta
ax.legend(loc='lower right',facecolor='white',frameon=True)
ax.set_xlabel('$m_1$')
ax.set_ylabel('$c_1$')
ax.set_title('value function in period 1')
ax.set_xlim([0,2])
ax.set_ylim([0,1.5]);
```

**Main takeaway:** The household lower its consumption when risk increases (such as in a recession). This is called **precautionary saving**.

## Simulation

**Step 1:** Solve and construct interpolators:


```python
Delta = 0.5
m1_vec,c1_vec,m2_vec,c2_vec = solve(par)

c1_interp = interpolate.RegularGridInterpolator([m1_vec], c1_vec,
                                                bounds_error=False,fill_value=None)

c2_interp = interpolate.RegularGridInterpolator([m2_vec], c2_vec,
                                                bounds_error=False,fill_value=None)
```


**Step 2:** Draw initail distribution of $m_1$ and simulate forward


```python
# a. draw initial m1
simN = 10000
sim_m1 = np.fmax(np.random.normal(1,0.1,size=simN),0) # "randomly" chosen distribution

# b. period 1
sim_c1 = c1_interp(sim_m1)
sim_a1 = sim_m1-sim_c1

# c. transition to period 2 with random draw
sim_m2 = (1+par.r)*sim_a1+np.random.choice([0.5,1.5],p=[0.5,0.5])

# d. period 2
sim_c2 = c2_interp(sim_m2)
```

**Step 3:** Plot distributions


```python
fig = plt.figure()
ax = fig.add_subplot(1,1,1)

ax.hist(sim_c1,bins=100,label='period 1')
ax.hist(sim_c2,bins=100,label='period 2')

ax.legend(loc='lower right',facecolor='white',frameon=True)
ax.set_xlabel('$c_t$')
ax.set_ylabel('freq.')
ax.set_title('consumption');
```

**Conclusion:** You can now solve models with complex heterogeneity and uncertainty, and simulate the implied dynamics. By introducing various policies you can quantify their effect not just for the average, but for the full distribution.

# Summary

**This lecture:**

1. Solving multidimensional optimization problems with and without gradients (and hessians)
2. Using multistart to alleviate problems with local minima (due to non-convexities)
3. Using penalty terms to solve constrained optimization problems 
4. Linear interpolation between known points
5. Solving dynamic optimization problems backwards period-by-period 

**Dynamic optimization:** Extremely useful technique. Can handle multiple periods, multiple states and choices, more shocks etc. You can solve general equilibrium models where the households solve such problems.

**Need more dynamic optimization?** [Mini-Course on Dynamic Programming](https://github.com/NumEconCopenhagen/ConsumptionSavingNotebooks/tree/master/00.%20DynamicProgramming#mini-course-in-dynamic-programming)

**Next lecture:** Canonical Economic Models.

2e3b00eadcb111532909703cda6dc0e86e0cc11c
ipynb
Jupyter Notebook
11/Numerical_optimization.ipynb
tarabrab/lectures-2021
41001cbce790c59e5c1f4ba19a744fdf9cbd952e
MIT
11/Numerical_optimization.ipynb
tarabrab/lectures-2021
41001cbce790c59e5c1f4ba19a744fdf9cbd952e
MIT
11/Numerical_optimization.ipynb
tarabrab/lectures-2021
41001cbce790c59e5c1f4ba19a744fdf9cbd952e
MIT
Qwen/Qwen-72B
1. YES
2. YES