---
title: "Optimization (cont.)"
subtitle: "Lecture 14"
author: "Dr. Colin Rundel"
footer: "Sta 663 - Spring 2025"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
    self-contained: true
execute: 
  echo: true
---

```{python setup}
#| include: false
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

import timeit

plt.rcParams['figure.dpi'] = 200
pd.set_option('display.width', 120)

from scipy import optimize
```

```{r r_setup}
#| include: false
knitr::opts_chunk$set(
  fig.align="center",
  cache=FALSE
)

local({
  hook_err_old <- knitr::knit_hooks$get("error")  # save the old hook
  knitr::knit_hooks$set(error = function(x, options) {
    # now do whatever you want to do with x, and pass
    # the new x to the old hook
    x = sub("## \n## Detailed traceback:\n.*$", "", x)
    x = sub("Error in py_call_impl\\(.*?\\)\\: ", "", x)
    hook_err_old(x, options)
  })
  
  hook_warn_old <- knitr::knit_hooks$get("warning")  # save the old hook
  knitr::knit_hooks$set(warning = function(x, options) {
    x = sub("<string>:1: ", "", x)
    hook_warn_old(x, options)
  })
})

options(width=90)
```

```{python utility}
#| include: false

# Code from https://scipy-lectures.org/ on optimization
def mk_quad(epsilon, ndim=2):
  def f(x):
    x = np.asarray(x)
    y = x.copy()
    y *= np.power(epsilon, np.arange(ndim))
    return .33*np.sum(y**2)
  
  def gradient(x):
    x = np.asarray(x)
    y = x.copy()
    scaling = np.power(epsilon, np.arange(ndim))
    y *= scaling
    return .33*2*scaling*y
  
  def hessian(x):
    scaling = np.power(epsilon, np.arange(ndim))
    return .33*2*np.diag(scaling)
  
  return f, gradient, hessian

def mk_rosenbrock(y=None):
  def f(x):
    x = np.asarray(x)
    y = 4*x
    y[0] += 1
    y[1:] += 3
    return np.sum(.5*(1 - y[:-1])**2 + (y[1:] - y[:-1]**2)**2)
  
  def gradient(x):
    x = np.asarray(x)
    y = 4*x
    y[0] += 1
    y[1:] += 3
    xm = y[1:-1]
    xm_m1 = y[:-2]
    xm_p1 = y[2:]
    der = np.zeros_like(y)
    der[1:-1] = 2*(xm - xm_m1**2) - 4*(xm_p1 - xm**2)*xm - .5*2*(1 - xm)
    der[0] = -4*y[0]*(y[1] - y[0]**2) - .5*2*(1 - y[0])
    der[-1] = 2*(y[-1] - y[-2]**2)
    return 4*der
  
  def hessian(x):
    x = np.asarray(x)
    y = 4*x
    y[0] += 1
    y[1:] += 3
    
    H = np.diag(-4*y[:-1], 1) - np.diag(4*y[:-1], -1)
    diagonal = np.zeros_like(y)
    diagonal[0] = 12*y[0]**2 - 4*y[1] + 2*.5
    diagonal[-1] = 2
    diagonal[1:-1] = 3 + 12*y[1:-1]**2 - 4*y[2:]*.5
    H = H + np.diag(diagonal)
    return 4*4*H
  
  return f, gradient, hessian

def super_fmt(value):
    if value > 1:
        if np.abs(int(value) - value) < .1:
            out = '$10^{%.1i}$' % value
        else:
            out = '$10^{%.1f}$' % value
    else:
        value = np.exp(value - .01)
        if value > .1:
            out = '%1.1f' % value
        elif value > .01:
            out = '%.2f' % value
        else:
            out = '%.2e' % value
    return out
```

```{python}
#| include: false
def plot_2d_traj(x, y, f, traj=None, title="", figsize=(5,5)):
  x_min, x_max = x
  y_min, y_max = y
  
  plt.figure(figsize=figsize, layout="constrained")
  
  x, y = np.mgrid[x_min:x_max:100j, y_min:y_max:100j]
  x = x.T
  y = y.T
  
  plt.figure(figsize=figsize)
  #plt.clf()
  #plt.axes([0, 0, 1, 1])
  
  X = np.concatenate((x[np.newaxis, ...], y[np.newaxis, ...]), axis=0)
  z = np.apply_along_axis(f, 0, X)
  log_z = np.log(z + .01)
  plt.imshow(
    log_z,
    extent=[x_min, x_max, y_min, y_max],
    cmap=plt.cm.gray_r, origin='lower',
    vmax=log_z.min() + 1.5*np.ptp(log_z)
  )
  contours = plt.contour(
    log_z,
    extent=[x_min, x_max, y_min, y_max],
    cmap=plt.cm.gnuplot, origin='lower'
  )
  
  plt.clabel(contours, inline=1, fmt=super_fmt, fontsize=12)
  
  if not traj is None:
    tx, ty = zip(*traj["x"])

    plt.plot(tx, ty, ".-b", ms = 10)
    plt.plot(tx[0], ty[0], ".r", ms = 15)
    plt.plot(tx[-1], ty[-1], ".c", ms = 15)
  
  if not title == "":
    plt.title(title)
  
  plt.xlim(x_min, x_max)
  plt.ylim(y_min, y_max)
  
  plt.show()
  plt.close('all')

```

```{python newtons}
#| include: false
def newtons_method(
  x, f, grad, hess, max_iter=100, max_back=10, tol=1e-8,
  alpha=0.5, beta=0.75
):
    res = {
      "x": [x], 
      "f": [f],
      "nit": 0,
      "nfev": 0,
      "njev": 0,
      "nhev": 0,
      "success": False
    }
    
    for i in range(max_iter):
      grad_f = grad(x)
      step = - np.linalg.solve(hess(x), grad_f) 

      res["njev"] += 1
      res["nhev"] += 1
      
      t = 1
      for j in range(max_back):
        res["nfev"] += 2
        if f(x+t*step) < f(x) + alpha * t * grad_f @ step:
          break
        t = t * beta
      
      x = x + t * step

      res["nit"] += 1

      if np.sqrt(np.sum((x - res["x"][-1])**2)) < tol:
        break
      
      res["x"].append(x)
      res["f"].append(f(x))
    
    return res
```

```{python cg}
#| include: false

def conjugate_gradient(x, f, grad, hess, max_iter=100, tol=1e-8):
    res = {
      "x": [x], 
      "f": [f],
      "nit": 0,
      "nfev": 0,
      "njev": 0,
      "nhev": 0,
      "success": False
    }

    r = grad(x)
    p = -r
    
    res["njev"] +=1

    for i in range(max_iter):
      H = hess(x)
      a = - r.T @ p / (p.T @ H @ p)
      x = x + a * p
      r = grad(x)
      b = (r.T @ H @ p) / (p.T @ H @ p)
      p = -r + b * p
      
      res["njev"] += 1
      res["nhev"] += 1
      res["nit"] += 1

      if np.sqrt(np.sum(r**2)) < tol:
        res["success"] = True
        break

      res["x"].append(x) 
      res["f"].append(f(x))
  
    return res
```

# Method Variants

## Method: CG in scipy

Scipy's optimize module implements the conjugate gradient algorithm using a variant proposed by Polak and Ribiere, this version does not use the Hessian when calculating the next step. The specific changes are:

* $\alpha_k$ is calculated via a line search along the direction $p_k$

* and the $\beta_{k+1}$ calculation is replaced as follows

$$
\beta_{k+1} = \frac{ r^T_{k+1} \, \nabla^2 f(x_k)  \, p_{k} }{p_k^T \, \nabla^2 f(x_k) \, p_k} 
\qquad \Rightarrow \qquad
\beta_{k+1}^{PR} = \frac{\nabla f(x_{k+1}) \left(\nabla f(x_{k+1}) - \nabla f(x_{k})\right)}{\nabla f(x_k)^T \, \nabla f(x_k)}
$$


## Method: Newton-CG & BFGS

These are both variants of Newton's method but do not require the Hessian (but can be used by the former if provided).

* Newton-Conjugate Gradient (Netwon-CG) algorithm uses a conjugate gradient algorithm to (approximately) invert the local Hessian

* The Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm iteratively approximates the inverse Hessian

    * Gradient is also not required and can similarly be approximated using finite differences
    


## Method: Nelder-Mead

This is a gradient free method that uses a series of simplexes which are used to iteratively bracket the minimum.

<iframe data-src="https://rundel.github.io/nelder-mead/" width="100%" height="450px" style="border:1px solid;border-radius: 5px;" data-external="1">
</iframe>

::: {.aside}
This js implementation comes from [greg-rychlewski/nelder-mead](https://github.com/greg-rychlewski/nelder-mead)
:::


## Method Summary

<br/>

::: {.small}

| SciPy Method     | Description                                                    | Gradient | Hessian  |
|:-----------------|:---------------------------------------------------------------|:--------:|:--------:|
|  ---             | Gradient Descent (naive w/ backtracking)                       |    ✓     |    ✗     |
|  ---             | Newton's method (naive w/ backtracking)                        |    ✓     |    ✓     |
|  ---             | Conjugate Gradient (naive)                                     |    ✓     |    ✓     |
| `"CG"`           | Nonlinear Conjugate Gradient (Polak and Ribiere variation)     |    ✓     |    ✗     |
| `"Newton-CG"`    | Truncated Newton method (Newton w/ CG step direction)          |    ✓     | Optional |
| `"BFGS"`         | Broyden, Fletcher, Goldfarb, and Shanno (Quasi-newton method)  | Optional |    ✗     |
| `"L-BFGS-B"`     | Limited-memory BFGS (Quasi-newton method)                      | Optional |    ✗     |
| `"Nelder-Mead"`  | Nelder-Mead simplex reflection method                          |    ✗     |    ✗     |
| 

:::

## Methods collection

::: {.small}
```{python}
def define_methods(x, f, grad, hess, tol=1e-8):
  return {
    "naive_newton":   lambda: newtons_method(x, f, grad, hess, tol=tol),
    "naive_cg":       lambda: conjugate_gradient(x, f, grad, hess, tol=tol),
    "CG":             lambda: optimize.minimize(f, x, jac=grad, method="CG", tol=tol),
    "newton-cg":      lambda: optimize.minimize(f, x, jac=grad, hess=None, method="Newton-CG", tol=tol),
    "newton-cg w/ H": lambda: optimize.minimize(f, x, jac=grad, hess=hess, method="Newton-CG", tol=tol),
    "bfgs":           lambda: optimize.minimize(f, x, jac=grad, method="BFGS", tol=tol),
    "bfgs w/o G":     lambda: optimize.minimize(f, x, method="BFGS", tol=tol),
    "l-bfgs-b":       lambda: optimize.minimize(f, x, method="L-BFGS-B", tol=tol),
    "nelder-mead":    lambda: optimize.minimize(f, x, method="Nelder-Mead", tol=tol)
  }
```
:::


## Method Timings

::: {.xsmall}
```{python}
x = (1.6, 1.1)
f, grad, hess = mk_quad(0.7)
methods = define_methods(x, f, grad, hess)
df = pd.DataFrame({
  key: timeit.Timer(methods[key]).repeat(10, 100) for key in methods
})
```
:::

```{python}
#| echo: false
g = sns.catplot(data=df.melt(), y="variable", x="value", aspect=1.8)
g.ax.set_xlabel("Time (sec)")
g.ax.set_ylabel("")
plt.show()
```

## Timings across cost functions

::: {.xsmall}
```{python}
def time_cost_func(n, x, name, cost_func, *args):
  f, grad, hess = cost_func(*args)
  methods = define_methods(x, f, grad, hess)
  
  return ( pd.DataFrame({
      key: timeit.Timer(
        methods[key]
      ).repeat(n, n) 
      for key in methods
    })
    .melt()
    .assign(cost_func = name)
  )

x = (1.6, 1.1)  

time_cost_df = pd.concat([
  time_cost_func(10, x, "Well-cond quad", mk_quad, 0.7),
  time_cost_func(10, x, "Ill-cond quad", mk_quad, 0.02),
  time_cost_func(10, x, "Rosenbrock", mk_rosenbrock)
])
```
:::


##

```{python}
#| echo: false
g = sns.catplot(data=time_cost_df, y="variable", x="value", hue="cost_func", alpha=0.5, aspect=2)
g.ax.set_xlabel("Time (sec)")
g.ax.set_ylabel("")
plt.show()
```


## Random starting locations

:::: {.tiny}
```{python}
pts = np.random.default_rng(seed=1234).uniform(-2,2, (20,2))

df = pd.concat([
  pd.concat([
    time_cost_func(3, x, "Well-cond quad", mk_quad, 0.7),
    time_cost_func(3, x, "Ill-cond quad", mk_quad, 0.02),
    time_cost_func(3, x, "Rosenbrock", mk_rosenbrock)
  ])
  for x in pts
])
```
::::

```{python}
#| echo: false
g = sns.catplot(data=df, y="variable", x="value", hue="cost_func", alpha=0.5, aspect=2)
g.ax.set_xlabel("Time (3 iter)")
g.ax.set_ylabel("")
plt.show()
``` 


## Profiling - BFGS (cProfile)

::: {.xsmall}
```{python}
import cProfile

f, grad, hess = mk_quad(0.7)
def run():
  for i in range(500):
    optimize.minimize(fun = f, x0 = (1.6, 1.1), jac=grad, method="BFGS", tol=1e-11)

cProfile.run('run()', sort="tottime")
```
:::

## Profiling - BFGS (pyinstrument)

::: {.xsmall}
```{python}
from pyinstrument import Profiler

f, grad, hess = mk_quad(0.7)

profiler = Profiler(interval=0.00001)

profiler.start()
opt = optimize.minimize(fun = f, x0 = (1.6, 1.1), jac=grad, method="BFGS", tol=1e-11)
p = profiler.stop()

profiler.write_html("Lec14_bfgs_quad.html")
```
:::

##

<iframe data-src="Lec14_bfgs_quad.html" width="100%" height="600px" style="border:1px solid;border-radius: 5px;" data-external="1">
</iframe>



## Profiling - Nelder-Mead

::: {.small}
```{python}
from pyinstrument import Profiler

f, grad, hess = mk_quad(0.7)

profiler = Profiler(interval=0.00001)

profiler.start()
opt = optimize.minimize(fun = f, x0 = (1.6, 1.1), method="Nelder-Mead", tol=1e-11)
p = profiler.stop()

profiler.write_html("Lec14_nm_quad.html")
```
:::


##

<iframe data-src="Lec14_nm_quad.html" width="100%" height="600px" style="border:1px solid;border-radius: 5px;" data-external="1">
</iframe>


## `optimize.minimize()` output

::: {.xsmall}
```{python}
f, grad, hess = mk_quad(0.7)
```
:::

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
optimize.minimize(
  fun = f, x0 = (1.6, 1.1), 
  jac=grad, method="BFGS"
)
```
:::

::: {.column width='50%'}
```{python}
optimize.minimize(
  fun = f, x0 = (1.6, 1.1), 
  jac=grad, hess=hess, method="Newton-CG"
)
```
:::
::::

## `optimize.minimize()` output (cont.)

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
optimize.minimize(
  fun = f, x0 = (1.6, 1.1), 
  jac=grad, method="CG"
)
```
:::

::: {.column width='50%'}
```{python}
optimize.minimize(
  fun = f, x0 = (1.6, 1.1), 
  jac=grad, method="Nelder-Mead"
)
```
:::
::::


## Collect

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
def run_collect(name, x0, cost_func, *args, tol=1e-8, skip=[]):
  f, grad, hess = cost_func(*args)
  methods = define_methods(x0, f, grad, hess, tol)
  
  res = []
  for method in methods:
    if method in skip: continue
    
    x = methods[method]()
    time = timeit.Timer(methods[method]).repeat(10, 25)

    d = {
      "name":    name,
      "method":  method,
      "nit":     x["nit"],
      "nfev":    x["nfev"],
      "njev":    x.get("njev"),
      "nhev":    x.get("nhev"),
      "success": x["success"],
      "time":    [time]
    }
    res.append( pd.DataFrame(d, index=[1]) )
  
  return pd.concat(res)
```
:::

::: {.column width='50%'}
```{python}
df = pd.concat([
  run_collect(
    name, (1.6, 1.1), 
    cost_func, 
    arg, 
    skip=['naive_newton', 'naive_cg']
  ) 
  for name, cost_func, arg in zip(
    ("Well-cond quad", "Ill-cond quad", "Rosenbrock"), 
    (mk_quad, mk_quad, mk_rosenbrock), 
    (0.7, 0.02, None)
  )
])
```
:::
::::

##

::: {.panel-tabset}

### Runtimes

```{python}
#| echo: false
g = sns.catplot(data=time_cost_df.query("variable != 'naive_newton' & variable != 'naive_cg'"), y="variable", x="value", hue="cost_func", alpha=0.5, aspect=2)
g.ax.set_xlabel("Time (10 iter)")
g.ax.set_ylabel("")
plt.show()
```

### Function calls

```{python}
#| echo: false
g = sns.catplot(
  y = "method", x = "value", hue = "variable", col="name", kind="bar",
  data = df.melt(id_vars=["name","method"], 
  value_vars=["nit", "nfev", "njev", "nhev"]).astype({"value": "float64"}),
  height=5, aspect=2./3.
).set_xlabels("n")
```

:::


## Exercise 1

Try minimizing the following function using different optimization methods starting from $x_0 = [0,0]$, which method(s) appear to work best?

$$
\begin{align}
f(x) = \exp(x_1-1) + \exp(-x_2+1) + (x_1-x_2)^2
\end{align}
$$
```{python}
#| echo: false
#| out-width: 40%
f = lambda x: np.exp(x[0]-1) + np.exp(-x[1]+1) + (x[0]-x[1])**2
plot_2d_traj((-2,3), (-2,3), f)
```

# MVN Example

## MVN density cost function

:::: {.columns}
::: {.column width='50%' .small}

For an $n$-dimensional multivariate normal we define <br/>
the $n \times 1$ vectors $x$ and $\mu$ and the $n \times n$ <br/>
covariance matrix $\Sigma$,

$$
\begin{align}
f(x) =& \det(2\pi\Sigma)^{-1/2} \\
      & \exp \left[-\frac{1}{2} (x-\mu)^T \Sigma^{-1} (x-\mu) \right] \\
\\
\nabla f(x) =& -f(x) \Sigma^{-1}(x-\mu) \\
\\
\nabla^2 f(x) =& f(x) \left( \Sigma^{-1}(x-\mu)(x-\mu)^T\Sigma^{-1} - \Sigma^{-1}\right) \\
\end{align}
$$

Our goal will be to find the mode (maximum) of this density.
:::

::: {.column width='50%' .xsmall}
```{python}
def mk_mvn(mu, Sigma):
  Sigma_inv = np.linalg.inv(Sigma)
  norm_const = 1 / (np.sqrt(np.linalg.det(2*np.pi*Sigma)))
  
  # Returns the negative density (since we want the max not min)
  def f(x):
    x_m = x - mu
    return -(norm_const * 
      np.exp( 
        -0.5 * x_m.T @ Sigma_inv @ x_m
      ) 
    ).item()
  
  def grad(x):
    return (-f(x) * Sigma_inv @ (x - mu))
  
  def hess(x):
    n = len(x)
    x_m = x - mu
    return f(x) * (
      (Sigma_inv @ x_m).reshape((n,1)) 
      @ (x_m.T @ Sigma_inv).reshape((1,n))
      - Sigma_inv
    )
  
  return f, grad, hess
```
:::
::::

::: {.aside}
From Section 8.1.1 of the [Matrix Cookbook](https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf)
:::


## Gradient checking

::: {.medium}
One of the most common issues when implementing an optimizer is to get the gradient calculation wrong which can produce problematic results. It is possible to numerically check the gradient function by comparing results between the gradient function and finite differences from the objective function via `optimize.check_grad()`.
:::


:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
# 5d
f, grad, hess = mk_mvn(
  np.zeros(5), np.eye(5,5)
)
```

```{python}
optimize.check_grad(f, grad, np.zeros(5))
optimize.check_grad(f, grad, np.ones(5))
```
:::

::: {.column width='50%'}
```{python}
# 10d
f, grad, hess = mk_mvn(
  np.zeros(10), np.eye(10)
)
```

```{python}
optimize.check_grad(f, grad, np.zeros(10))
optimize.check_grad(f, grad, np.ones(10))
```
:::
::::

## Gradient checking (wrong gradient)

::: {.xsmall}
```{python}
wrong_grad = lambda x: 5*grad(x)
```
:::

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
# 5d
f, grad, hess = mk_mvn(
  np.zeros(5), np.eye(5)
)
```

```{python}
optimize.check_grad(f, wrong_grad, np.zeros(5))
optimize.check_grad(f, wrong_grad, np.ones(5))
```
:::

::: {.column width='50%'}
```{python}
# 10d
f, grad, hess = mk_mvn(
  np.zeros(10), np.eye(10)
)
```

```{python}
optimize.check_grad(f, wrong_grad, np.zeros(10))
optimize.check_grad(f, wrong_grad, np.ones(10))
```
:::
::::

. . .

<br/>

::: {.center .large}
Why does `np.ones()` detect an issue but `np.zeros()` does not?
:::


## Hessian checking

Note since the gradient of the gradient is the hessian we can use this function to check our implementation of the hessian as well, just use `grad()` as `func` and `hess()` as `grad`.

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
# 5d
f, grad, hess = mk_mvn(np.zeros(5), np.eye(5))
```

```{python}
optimize.check_grad(grad, hess, np.zeros(5))
optimize.check_grad(grad, hess, np.ones(5))
```
:::

::: {.column width='50%'}
```{python}
# 10d
f, grad, hess = mk_mvn(np.zeros(10), np.eye(10))
```

```{python}
optimize.check_grad(grad, hess, np.zeros(10))
optimize.check_grad(grad, hess, np.ones(10))
```
:::
::::


## Unit MVNs

::: {.xsmall}
```{python}
rng = np.random.default_rng(seed=1234)
runif = rng.uniform(-1,1, size=25)

df = pd.concat([
  run_collect(
    name, runif[:n], mk_mvn, 
    np.zeros(n), np.eye(n), 
    tol=1e-10, 
    skip=['naive_newton', 'naive_cg', 'bfgs w/o G', 'newton-cg w/ H', 'l-bfgs-b', 'nelder-mead']
  ) 
  for name, n in zip(
    ("5d", "10d", "25d"), 
    (5, 10, 25)
  )
])
```
:::

## Performance (Unit MVNs)

::: {.panel-tabset}

### Run times

```{python}
#| echo: false
g = sns.catplot(
    data=(df[["name","method","time"]]
        .explode("time")
    ),
    y="method", 
    x="time", 
    hue="name", 
    alpha=0.5, 
    aspect=2
)
g.ax.set_xlabel("Time (secs)")
g.ax.set_ylabel("")
plt.show()
```

### Function calls

```{python}
#| echo: false
g = sns.catplot(
  y = "method", x = "value", hue = "variable", col="name", kind="bar",
  col_wrap=3,
  data = df.melt(
    id_vars=["name","method"], 
    value_vars=["nit", "nfev", "njev", "nhev"]).astype({"value": "float64"}
  ),
  aspect=0.75
).set_xlabels("n").set(xscale="log")
```

:::



## What's going on? (good)

::: {.xsmall}
```{python}
n = 5
f, grad, hess = mk_mvn(np.zeros(n), np.eye(n))
```
:::


:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
optimize.minimize(
  f, runif[:n], jac=grad, 
  method="newton-cg", tol=1e-10
)
```
:::

::: {.column width='50%'}
```{python}
optimize.minimize(
  f, runif[:n], jac=grad, 
  method="bfgs", tol=1e-10
)
```
:::
::::


## What's going on? (okay)

::: {.xsmall}
```{python}
n = 10
f, grad, hess = mk_mvn(np.zeros(n), np.eye(n))
```
:::


:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
optimize.minimize(
  f, runif[:n], jac=grad, 
  method="newton-cg", tol=1e-10
)
```
:::

::: {.column width='50%'}
```{python}
optimize.minimize(
  f, runif[:n], jac=grad, 
  method="bfgs", tol=1e-10
)
```
:::
::::


## What's going on? (bad)

::: {.xsmall}
```{python}
n = 25
f, grad, hess = mk_mvn(np.zeros(n), np.eye(n))
```
:::


:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
optimize.minimize(
  f, runif[:n], jac=grad, 
  method="newton-cg", tol=1e-10
)
```
:::

::: {.column width='50%'}
```{python}
optimize.minimize(
  f, runif[:n], jac=grad, 
  method="bfgs", tol=1e-10
)
```
:::
::::


## All bad?

::: {.xsmall}
```{python}
optimize.minimize(
  f, runif[:n], jac=grad, 
  method="nelder-mead", tol=1e-10
)
```
:::



## Options (newton-cg)

::: {.xsmall}
```{python}
optimize.show_options(solver="minimize", method="newton-cg")
```
:::

## Options (bfgs)

::: {.xsmall}
```{python}
optimize.show_options(solver="minimize", method="bfgs")
```
:::

## Options (Nelder-Mead)

::: {.xsmall}
```{python}
optimize.show_options(solver="minimize", method="newton-cg")
```
:::


## SciPy implementation

The following code comes from SciPy's `minimize()` implementation:

::: {.small}
```{python}
#| eval: false
if tol is not None:
  options = dict(options)
  if meth == 'nelder-mead':
      options.setdefault('xatol', tol)
      options.setdefault('fatol', tol)
  if meth in ('newton-cg', 'powell', 'tnc'):
      options.setdefault('xtol', tol)
  if meth in ('powell', 'l-bfgs-b', 'tnc', 'slsqp'):
      options.setdefault('ftol', tol)
  if meth in ('bfgs', 'cg', 'l-bfgs-b', 'tnc', 'dogleg',
              'trust-ncg', 'trust-exact', 'trust-krylov'):
      options.setdefault('gtol', tol)
  if meth in ('cobyla', '_custom'):
      options.setdefault('tol', tol)
  if meth == 'trust-constr':
      options.setdefault('xtol', tol)
      options.setdefault('gtol', tol)
      options.setdefault('barrier_tol', tol)
```
:::

::: {.aside}
See code in context on GitHub [here](https://github.com/scipy/scipy/blob/a438ba6ef4061c28d79657b525ed2378154dfea8/scipy/optimize/_minimize.py#L610-L627)
:::



## Fixing our tolerances?

::: {.xsmall}
```{python}
n = 25
f, grad, hess = mk_mvn(np.zeros(n), np.eye(n))
```
:::


:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
optimize.minimize(
  f, runif[:n], jac=grad, 
  method="newton-cg", tol=1e-16
)
```
:::

::: {.column width='50%'}
```{python}
optimize.minimize(
  f, runif[:n], jac=grad, 
  method="bfgs", tol=1e-16
)
```
:::
::::


## Limits of floating point precision

Every type of floating point value has finite precision due to the limitations of how it is represented. This value is typically refered to as the machine epsilon value, this is the smallest possible spacing between 1.0 and the next representable floating-point number.

:::: {.columns .small}
::: {.column width='50%'}
```{python}
np.finfo(np.float64).eps
np.finfo(np.float32).eps
np.finfo(np.float16).eps
```

:::

::: {.column width='50%'}
```{python}
1+np.finfo(np.float64).eps > 1
1+np.finfo(np.float64).eps/2 > 1
```
:::
::::


## Fixes?

::: {.xsmall}
```{python}
def mk_prop_mvn(mu, Sigma):
  Sigma_inv = np.linalg.inv(Sigma)
  #norm_const = 1 / (np.sqrt(np.linalg.det(2*np.pi*Sigma)))
  norm_const = 1
  
  # Returns the negative density (since we want the max not min)
  def f(x):
    x_m = x - mu
    return -(norm_const * 
      np.exp( 
        -0.5 * x_m.T @ Sigma_inv @ x_m
      ) 
    ).item()
  
  def grad(x):
    return (-f(x) * Sigma_inv @ (x - mu))
  
  def hess(x):
    n = len(x)
    x_m = x - mu
    return f(x) * (
      (Sigma_inv @ x_m).reshape((n,1)) 
      @ (x_m.T @ Sigma_inv).reshape((1,n))
      - Sigma_inv
    )
  
  return f, grad, hess
```
:::

## 

::: {.xsmall}
```{python}
n = 25
f, grad, hess = mk_prop_mvn(np.zeros(n), np.eye(n))
```
:::


:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
optimize.minimize(
  f, runif[:n], jac=grad, 
  method="newton-cg", tol=1e-10
)
```
:::

::: {.column width='50%'}
```{python}
optimize.minimize(
  f, runif[:n], jac=grad, 
  method="bfgs", tol=1e-10
)
```
:::
::::


## Performance 

```{python}
#| include: false
rng = np.random.default_rng(seed=1234)
runif = rng.uniform(-1,1, size=25)

df = pd.concat([
  run_collect(
    name, runif[:n], mk_prop_mvn, 
    np.zeros(n), np.eye(n), 
    tol=1e-10, 
    skip=['naive_newton', 'naive_cg', 'bfgs w/o G', 'newton-cg w/ H', 'l-bfgs-b', 'nelder-mead']
  ) 
  for name, n in zip(
    ("5d", "10d", "25d"), 
    (5, 10, 25)
  )
])
```


::: {.panel-tabset}

### Run times

```{python}
#| echo: false
g = sns.catplot(
    data=(df[["name","method","time"]]
        .explode("time")
    ),
    y="method", 
    x="time", 
    hue="name", 
    alpha=0.5, 
    aspect=2
)
g.ax.set_xlabel("Time (secs)")
g.ax.set_ylabel("")
plt.show()
```

### Function calls

```{python}
#| echo: false
g = sns.catplot(
  y = "method", x = "value", hue = "variable", col="name", kind="bar",
  col_wrap=3,
  data = df.melt(
    id_vars=["name","method"], 
    value_vars=["nit", "nfev", "njev", "nhev"]).astype({"value": "float64"}
  ),
  aspect=0.75
).set_xlabels("n").set(xscale="log")
```

:::






## Some general advice

* Having access to the gradient is almost always helpful / necessary

* Having access to the hessian can be helpful, but usually does not significantly improve things

* The curse of dimensionality is real 

  * Be careful with `tol` - it means different things for different methods
  
* In general, **BFGS** or **L-BFGS** should be a first choice for most problems (either well- or ill-conditioned)

  * **CG** can perform better for well-conditioned problems with cheap function evaluations


# Maximum Likelihood example

## Normal MLE

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
from scipy.stats import norm

n = norm(-3.2, 1.25)
x = n.rvs(size=100, random_state=1234)
{'µ': x.mean(), 'σ': x.std()}
```
```{python}
mle_norm = lambda θ: -np.sum(
  norm.logpdf(x, loc=θ[0], scale=θ[1])
)
```
:::

::: {.column width='50%' .fragment}
```{python}
mle_norm([0,1])
mle_norm([-3, 1])
mle_norm([-3.2, 1.25])
mle_norm([-3.3, 1.25])
```
:::
::::


## Minimizing

::: {.small}
```{python}
optimize.minimize(mle_norm, x0=[0,1], method="bfgs")
```
:::

## Adding constraints

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
def mle_norm2(θ): 
  if θ[1] <= 0:
    return np.inf
  else:
    return -np.sum(
      norm.logpdf(x, loc=θ[0], scale=θ[1])
    )
```
:::

::: {.column width='50%' .fragment}
```{python}
optimize.minimize(mle_norm2, x0=[0,1], method="bfgs")
```
:::
::::

. . .

::: {.xsmall}
```{python}
{'µ': x.mean(), 'σ': x.std()}
```
:::


## Specifying Bounds

It is also possible to specify bounds via `bounds` but this is only available for certain optimization methods (i.e. Nelder-Mead & L-BFGS-B).

::: {.small}
```{python}
optimize.minimize(
  mle_norm, x0=[0,1], method="l-bfgs-b",
  bounds = [(-1e16, 1e16), (1e-16, 1e16)]
)
```
:::


## Exercise 2

Using `optimize.minimize()` recover the shape and scale parameters for these data using MLE.

::: {.small}
```{python}
from scipy.stats import gamma

g = gamma(a=2.0, scale=2.0)
x = g.rvs(size=100, random_state=1234)
x.round(2)
```
:::
