---
title: "Numerical optimization"
subtitle: "Lecture 12"
author: "Dr. Colin Rundel"
footer: "Sta 663 - Spring 2023"
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

plt.rcParams['figure.dpi'] = 200

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)
  })
})
```

## Numerical optimization - line search

Today we will be discussing one particular approach for numerical optimization - line search. This is a family of algorithmic approaches that attempt to find (global or local) minima via iteration on an initial guess. Generally they are an attempt to solve,

$$
\underset{\alpha>0}{\text{min}} f(x_k + \alpha \, p_k)
$$
where $f()$ is the function we are attempting to minimize, $x_k$ is our current guess at iteration $k$ and $\alpha$ is the step length and $p_k$ is the direction of movement.


We will only be dipping our toes in the water of this area but the goal is to provide some context for some of the more common (and easier) use cases.


## Naive Gradient Descent

We will start with a naive approach to gradient descent where we choose a fixed step size and determine the direction based on the gradient of the function at each iteration.


::: {.small}
```{python}
def grad_desc_1d(x0, f, grad, step, max_step=100, tol = 1e-6):
  all_x_i = [x0]
  all_f_i = [f(x0)]
  x_i = x0
  
  try:
    for i in range(max_step):
      dx_i = grad(x_i)
      x_i = x_i - dx_i * step
      f_x_i = f(x_i)
      all_x_i.append(x_i)
      all_f_i.append(f_x_i)
      
      if np.abs(dx_i) < tol: break
    
  except OverflowError as err:
    print(f"{type(err).__name__}: {err}")
  
  if len(all_x_i) == max_step+1:
    print("Warning - Failed to converge!")
  
  return all_x_i, all_f_i
```
:::

```{python include=FALSE}
def plot_1d_traj(x, f, traj, title="", figsize=(5,3)):
  plt.figure(figsize=figsize, layout="constrained")
  
  x_range = x[1]-x[0]

  x_focus = np.linspace(x[0], x[1], 101)
  x_ext = np.linspace(x[0]-0.2*x_range, x[1]+0.2*x_range, 141)

  plt.plot(x_focus, f(x_focus), "-k")
  
  xlim = plt.xlim()
  ylim = plt.ylim()
  
  plt.plot(x_ext, f(x_ext), "-k")

  plt.plot(traj[0], traj[1], ".-b", ms = 10)
  plt.plot(traj[0][0], traj[1][0], ".r", ms = 15)
  plt.plot(traj[0][-1], traj[1][-1], ".c", ms = 15)

  plt.xlim(xlim)
  plt.ylim(ylim)
  
  plt.show()
  
  plt.close('all')
```

## A basic example

:::: {.columns .small}
::: {.column width='50%'}
$$
\begin{aligned}
f(x) &= x^2 \\
\nabla f(x) &= 2x
\end{aligned}
$$
:::

::: {.column width='50%'}
```{python}
f = lambda x: x**2
grad = lambda x: 2*x
```
:::
::::

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 90%
opt = grad_desc_1d(-2., f, grad, step=0.25)
plot_1d_traj( (-2, 2), f, opt )
```
:::

::: {.column width='50%' .fragment}
```{python}
#| out-width: 90%
opt = grad_desc_1d(-2., f, grad, step=0.5)
plot_1d_traj( (-2, 2), f, opt )
```
:::
::::


## Where can it go wrong?

If you pick a bad step size then bad things can happen,

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 90%
opt = grad_desc_1d(-2, f, grad, step=0.9)
plot_1d_traj( (-2,2), f, opt )
```
:::

::: {.column width='50%' .fragment}
```{python}
#| out-width: 90%
opt = grad_desc_1d(-2, f, grad, step=1)
plot_1d_traj( (-2,2), f, opt )
```
:::
::::



## Local minima of a quartic

Since the function is no longer convex - both starting point and step size matter.

:::: {.columns .small}
::: {.column width='50%'}
$$
\begin{aligned}
f(x) &= x^4 + x^3 -x^2 - x \\
\nabla f(x) &= 4x^3 + 3x^2 - 2x - 1
\end{aligned}
$$
:::

::: {.column width='50%'}
```{python}
f = lambda x: x**4 + x**3 - x**2 - x 
grad = lambda x: 4*x**3 + 3*x**2 - 2*x - 1
```
:::
::::

. . .

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

```{python out.width="90%", error=TRUE}
opt = grad_desc_1d(-1.5, f, grad, step=0.2)
plot_1d_traj( (-1.5, 1.5), f, opt )
```
:::

::: {.column width='50%' .fragment}
```{python out.width="90%", error=TRUE}
opt = grad_desc_1d(-1.5, f, grad, step=0.25)
plot_1d_traj( (-1.5, 1.5), f, opt)
```
:::
::::


## Alternative starting points

:::: {.columns .small}
::: {.column width='50%'}
```{python out.width="90%", error=TRUE}
opt = grad_desc_1d(1.5, f, grad, step=0.2)
plot_1d_traj( (-1.75, 1.5), f, opt )
```
:::

::: {.column width='50%' .fragment}
```{python out.width="90%", error=TRUE}
opt = grad_desc_1d(1.25, f, grad, step=0.2)
plot_1d_traj( (-1.75, 1.5), f, opt)
```
:::
::::


## Problematic step sizes

If the step size is too large it is possible for the algorithm to 

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

```{python out.width="90%", error=TRUE}
opt = grad_desc_1d(-1.5, f, grad, step=0.75)
plot_1d_traj( (-2, 2), f, opt )
```
:::

::: {.column width='50%' .fragment}
```{python out.width="90%", error=TRUE}
opt = grad_desc_1d(1.5, f, grad, step=0.25)
plot_1d_traj( (-2, 2), f, opt)
```
:::
::::

##  Gradient Descent w/ backtracking

:::: {.columns .small}
::: {.column width='50%'}
As we have just seen having too large of a step can 
be problematic, one solution is to allow the step size
to adapt.

Backtracking involves checking if the proposed move is
advantageous (i.e. $f(x_k+\alpha p_k) < f(x_k)$),

* If it is advantageous then accept
  $x_{k+1} = x_k+\alpha p_k$.

* If not, shrink $\alpha$ by a factor $\tau$ (e.g. 0.5)
  and check again.
  
Pick larger $\alpha$ to start as this will not fix
inefficiency of small step size.
:::

::: {.column width='50%' .fragment}
```{python}
#| code-line-numbers: 1,9-14
def grad_desc_1d_bt(
  x, f, grad, step, tau=0.5, 
  max_step=100, max_back=10, tol = 1e-6
):
  all_x_i = [x]
  all_f_i = [f(x)]
  
  try:
    for i in range(max_step):
      dx = grad(x)
      
      for j in range(max_back):
        new_x = x + step * (-dx)
        new_f_x = f(new_x)
        if (new_f_x < all_f_i[-1]): 
          break
        step = step * tau
      
      x = new_x
      f_x = new_f_x
      all_x_i.append(x)
      all_f_i.append(f_x)
      
      if np.abs(dx) < tol:
        break
    
  except OverflowError as err:
    print(f"{type(err).__name__}: {err}")
  
  if len(all_x_i) == max_step+1:
    print("Warning - Failed to converge!")
  
  return all_x_i, all_f_i
```
:::
::::

::: {.aside}
This is a simplified (hand wavy) version of the [Armijo-Goldstein condition](https://en.wikipedia.org/wiki/Backtracking_line_search) <br/>
Check $f(x_k-\alpha \nabla f(x_k)) \leq f(x_k) - c \alpha (\nabla f(x_k))^2$.
:::

##

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| error: true
#| out-width: 90%
opt = grad_desc_1d_bt(-1.5, f, grad, 
                      step=0.75, tau=0.5)
plot_1d_traj( (-1.5, 1.5), f, opt )
```
:::

::: {.column width='50%' .fragment}
```{python}
#| error: true
#| out-width: 90%
opt = grad_desc_1d_bt(1.5, f, grad, 
                      step=0.25, tau=0.5)
plot_1d_traj( (-1.5, 1.5), f, opt)
```
:::
::::




## A 2d cost function

```{python 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():
  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
```


We will be using `mk_quad()` to create quadratic functions with varying conditioning (as specified by the `epsilon` parameter).

$$
\begin{align}
f(x,y) &= 0.33(x^2 + \epsilon^2 y^2 ) \\
\nabla f(x,y) &= \left[ \begin{matrix}
0.66 \, x \\
0.66 \, \epsilon^2 \, y 
\end{matrix} \right] \\
\nabla^2 f(x,y) &= \left[\begin{array}{cc}
0.66  & 0 \\
0     & 0.66 \, \epsilon^2
\end{array}\right]
\end{align}
$$


## Examples

```{python include=FALSE}
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

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*log_z.ptp()
  )
  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:
    plt.plot(traj[0], traj[1], ".-b", ms = 10)
    plt.plot(traj[0][0], traj[1][0], ".r", ms = 15)
    plt.plot(traj[0][-1], traj[1][-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')
  
```

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 100%
f, grad, hess = mk_quad(0.7)
plot_2d_traj((-1,2), (-1,2), f, 
             title="$\\epsilon=0.7$")
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 100%
f, grad, hess = mk_quad(0.05)
plot_2d_traj((-1,2), (-1,2), f, 
             title="$\\epsilon=0.05$")
```
:::
::::

## 2d gradient descent w/ backtracking

::: {.small}
```{python}
def grad_desc_2d(x0, f, grad, step, tau=0.5, max_step=100, max_back=10, tol = 1e-6):
  x_i = x0
  all_x_i = [x_i[0]]
  all_y_i = [x_i[1]]
  all_f_i = [f(x_i)]
  
  for i in range(max_step):
    dx_i = grad(x_i)
    
    for j in range(max_back):
      new_x_i = x_i - dx_i * step
      new_f_i = f(new_x_i)
      if (new_f_i < all_f_i[-1]): break
      step = step * tau
      
    x_i, f_i = new_x_i, new_f_i
    all_x_i.append(x_i[0])
    all_y_i.append(x_i[1])
    all_f_i.append(f_i)
      
    if np.sqrt(np.sum(dx_i**2)) < tol:
      break
  
  return all_x_i, all_y_i, all_f_i
```
:::


## Well conditioned cost function

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.7)
opt = grad_desc_2d((1.6, 1.1), f, grad, step=1)
plot_2d_traj((-1,2), (-1,2), f, 
             title="$\\epsilon=0.7$", traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.7)
opt = grad_desc_2d((1.6, 1.1), f, grad, step=2)
plot_2d_traj((-1,2), (-1,2), f, 
             title="$\\epsilon=0.7$", traj=opt)
```
:::
::::


## Ill-conditioned cost function

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.05)
opt = grad_desc_2d((1.6, 1.1), f, grad, step=1)
plot_2d_traj((-1,2), (-1,2), f, 
             title="$\\epsilon=0.05$", traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.05)
opt = grad_desc_2d((1.6, 1.1), f, grad, step=2)
plot_2d_traj((-1,2), (-1,2), f, 
             title="$\\epsilon=0.05$", traj=opt)
```
:::
::::


## Rosenbrock function (very ill conditioned)

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = grad_desc_2d((1.6, 1.1), f, grad, step=0.25)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = grad_desc_2d((-0.5, 0), f, grad, step=0.25)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::
::::


## Taylor Expansion

For any arbitrary smooth function, we can construct a 2nd order Taylor approximation as follows,

$$
\begin{align}
f(x_k + \alpha \, p_k) 
&= f(x_k) + \alpha \, p_k^T \nabla f(x_k + \alpha \, p_k) \\
&= f(x_k) + \alpha \, p_k^T \nabla f(x_k) + \frac{1}{2} \alpha^2 p_k^T \, \nabla^2 f(x_k + \alpha \, p_k) \, p_k \\
&\approx f(x_k) + \alpha \, p_k^T \nabla f(x_k) + \frac{1}{2} \alpha^2 p_k^T \, \nabla^2 f(x_k) \, p_k
\end{align}
$$


## Newton's Method in 1d

Lets simplify things for now and consider just the 1d case and write $\alpha\,p_k$ as $\Delta$,

$$
f(x_k + \Delta) \approx f(x_k) + \Delta f'(x_k) + \frac{1}{2} \Delta^2 f''(x_k)
$$

to find the $\Delta$ that minimizes this function we can take a derivative with regard to $\Delta$ and set the equation equal to zero which gives,

$$
0 = f'(x_k) + \Delta f''(x_k) \;\; \Rightarrow \;\; \Delta = -\frac{f'(x_k)}{f''(x_k)}
$$
which then suggests an iterative update rule of 

$$
x_{k+1} = x_{k} -\frac{f'(x_k)}{f''(x_k)}
$$


## Generalizing to $n$d

Based on the same argument we can see the follow result for a function in $\mathbb{R}^n$,

$$
f(x_k + \Delta) \approx f(x_k) + \Delta^T \nabla f(x_k) + \frac{1}{2} \Delta^T \, \nabla^2 f(x_k) \,\Delta
$$

$$
0 = \nabla f(x_k) + \nabla^2 f(x_k) \, \Delta \;\; \Rightarrow \;\; \Delta = -\left(\nabla^2 f(x_k)\right)^{-1} \nabla f(x_k) f(x_k)
$$
which then suggests an iterative update rule of 

$$
x_{k+1} = x_{k} - (\nabla^2 f(x_k))^{-1} \, \nabla f(x_k)
$$


##

::: {.small}
```{python}
#| code-line-numbers: 1,9-19
def newtons_method(x0, f, grad, hess, max_iter=100, max_back=10, tol=1e-8):
    all_x_i = [x0[0]]
    all_y_i = [x0[1]]
    all_f_i = [f(x0)]
    
    x_i = x0
    
    for i in range(max_iter):
      g_i = grad(x_i)
      step = - np.linalg.solve(hess(x_i), g_i)
      
      for j in range(max_back):
        new_x_i = x_i + step
        new_f_i = f(new_x_i)
      
        if (new_f_i < all_f_i[-1]):
          break
      
        step /= 2
      
      x_i, f_i = new_x_i, new_f_i
      
      all_x_i.append(x_i[0])
      all_y_i.append(x_i[1])
      all_f_i.append(f_i)
      
      if np.sqrt(np.sum(g_i**2)) < tol:
        break
    
    return all_x_i, all_y_i, all_f_i
```
:::


::: {.aside}
Based on Chapter 5.1 from [Core Statistics](https://www.maths.ed.ac.uk/~swood34/core-statistics.pdf)
:::


## Well conditioned quadratic cost function

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.7)
opt = newtons_method((1.6, 1.1), f, grad, hess)
plot_2d_traj((-1,2), (-1,2), f, traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.05)
opt = newtons_method((1.6, 1.1), f, grad, hess)
plot_2d_traj((-1,2), (-1,2), f, traj=opt)
```
:::
::::


## Rosenbrock function (very ill conditioned)

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = newtons_method((1.6, 1.1), f, grad, hess)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = newtons_method((-0.5, 0), f, grad, hess)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::
::::


## Conjugate gradients

This is a general approach for solving a system of linear equations with the form $Ax=b$ where $A$ is an $n \times n$ symmetric positive definite matrix and b is $n \times 1$ with $x$ unknown.

This type of problem can also be expressed as a quadratic minimization problems of the form,

$$
\underset{x}{\text{min}} \; f(x) = \frac{1}{2} x^T \, A \, x - b^T x + c
$$

The goal is then to find $n$ conjugate vectors ( $p^T_i \, A \, p_j = 0$ for all $i \neq j$) and their coefficients such that 

$$ x_* = \sum_{i=1}^n \alpha_i \, p_i $$


## Conjugate gradient algorithm

:::: {.columns .small}
::: {.column width='50%'}
Given $x_0$ we set the following initial values,

$$\begin{align}
r_0 &= \nabla f(x_0) \\
p_0 &= -r_0 \\
k &= 0
\end{align}$$

while $\|r_k\|_2 > \text{tol}$,

$$
\begin{align}
\alpha_k &= \frac{r_k^T \, p_k}{p_k^T \, \nabla^2 f(x_k) \, p_k} \\
x_{k+1} &= x_k + \alpha_k \, p_k \\
r_{k+1} &= \nabla f(x_{k+1}) \\
\beta_{k} &= \frac{ r^T_{k+1} \, \nabla^2 f(x_k)  \, p_{k} }{p_k^T \, \nabla^2 f(x_k) \, p_k} \\
p_{k+1} &= -r_{k+1} + \beta_{k} \, p_k \\
k &= k+1
\end{align}
$$
:::

::: {.column width='50%'}
```{python}
#| code-line-numbers: 1-2,11-18
def conjugate_gradient(x0, f, grad, hess, 
                       max_iter=100, tol=1e-8):
    all_x_i = [x0[0]]
    all_y_i = [x0[1]]
    all_f_i = [f(x0)]
    
    x_i = x0
    r_i = grad(x0)
    p_i = -r_i
    
    for i in range(max_iter):
      a_i = - r_i.T @ p_i / (p_i.T @ hess(x_i) @ p_i)
      x_i_new = x_i + a_i * p_i
      r_i_new = grad(x_i_new)
      b_i = (r_i_new.T @ hess(x_i) @ p_i) / (p_i.T @ hess(x_i) @ p_i)
      p_i_new = -r_i_new + b_i * p_i
      
      x_i, r_i, p_i = x_i_new, r_i_new, p_i_new
      
      all_x_i.append(x_i[0])
      all_y_i.append(x_i[1])
      all_f_i.append(f(x_i))
      
      if np.sqrt(np.sum(r_i_new**2)) < tol:
        break
    
    return all_x_i, all_y_i, all_f_i
```
:::
::::

::: {.aside}
From Chapter 5.1 of [Numerical Optimization](https://find.library.duke.edu/catalog/DUKE004973775) 2006
:::


## Trajectory

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.7)
opt = conjugate_gradient((1.6, 1.1), f, grad, hess)
plot_2d_traj((-1,2), (-1,2), f, title="$\\epsilon=0.7$", traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.05)
opt = conjugate_gradient((1.6, 1.1), f, grad, hess)
plot_2d_traj((-1,2), (-1,2), f, title="$\\epsilon=0.05$", traj=opt)
```
:::
::::


## Rosenbrock's function 

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = conjugate_gradient((1.6, 1.1), f, grad, hess)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = conjugate_gradient((-0.5, 0), f, grad, hess)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::
::::


## CG in scipy

Scipy's optimize module implements the conjugate gradient algorithm by Polak and Ribiere, a variant that does not require the Hessian,

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

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

* $\beta_{k+1}$ is replaced with

$$
\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)}
$$

:::

::: {.column width='50%'}
```{python}
#| code-line-numbers: 1,6-15
def conjugate_gradient_scipy(x0, f, grad, tol=1e-8):
    all_x_i = [x0[0]]
    all_y_i = [x0[1]]
    all_f_i = [f(x0)]
    
    def store(X):
        x, y = X
        all_x_i.append(x)
        all_y_i.append(y)
        all_f_i.append(f(X))
    
    optimize.minimize(
      f, x0, jac=grad, method="CG", 
      callback=store, tol=tol
    )
    
    return all_x_i, all_y_i, all_f_i
```
:::
::::


## Trajectory

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.7)
opt = conjugate_gradient_scipy((1.6, 1.1), f, grad)
plot_2d_traj((-1,2), (-1,2), f, title="$\\epsilon=0.7$", traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.05)
opt = conjugate_gradient_scipy((1.6, 1.1), f, grad)
plot_2d_traj((-1,2), (-1,2), f, title="$\\epsilon=0.05$", traj=opt)
```
:::
::::


## Rosenbrock's function 

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = conjugate_gradient_scipy((1.6, 1.1), f, grad)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = conjugate_gradient_scipy((-0.5, 0), f, grad)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::
::::


## Method: Newton-CG

Is a variant of Newtons method but does not require inverting the Hessian, or even a Hessian function (for the latter case it is estimated by finite differencing of the gradient)

:::: {.columns .small}
::: {.column width='50%'}
```{python}
def newton_cg(x0, f, grad, hess=None, tol=1e-8):
    all_x_i = [x0[0]]
    all_y_i = [x0[1]]
    all_f_i = [f(x0)]
    
    def store(X):
        x, y = X
        all_x_i.append(x)
        all_y_i.append(y)
        all_f_i.append(f(X))
    
    optimize.minimize(
      f, x0, jac=grad, hess=hess, tol=tol,
      method="Newton-CG", callback=store 
    )
    
    return all_x_i, all_y_i, all_f_i
```
:::
::::


## Trajectory - well conditioned

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.7)
opt = newton_cg((1.6, 1.1), f, grad)
plot_2d_traj((-1,2), (-1,2), f, traj=opt, title="w/o hessian")
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.7)
opt = newton_cg((1.6, 1.1), f, grad, hess)
plot_2d_traj((-1,2), (-1,2), f, traj=opt, title="w/ hessian")
```
:::
::::


## Trajectory - ill-conditioned

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.05)
opt = newton_cg((1.6, 1.1), f, grad)
plot_2d_traj((-1,2), (-1,2), f, traj=opt, title="w/o hessian")
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.05)
opt = newton_cg((1.6, 1.1), f, grad, hess)
plot_2d_traj((-1,2), (-1,2), f, traj=opt, title="w/ hessian")
```
:::
::::


## Rosenbrock's function 

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = newton_cg((1.6, 1.1), f, grad)
plot_2d_traj((-2,2), (-2,2), f, traj=opt, title="w/o hessian")
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = newton_cg((1.6, 1.1), f, grad, hess)
plot_2d_traj((-2,2), (-2,2), f, traj=opt, title="w/ hessian")
```
:::
::::


## Method: BFGS

The Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm is a quasi-newton which iterative improves its approximation of the Hessian,

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| code-line-numbers: 1,12-15
def bfgs(x0, f, grad, hess=None, tol=1e-8):
    all_x_i = [x0[0]]
    all_y_i = [x0[1]]
    all_f_i = [f(x0)]
    
    def store(X):
        x, y = X
        all_x_i.append(x)
        all_y_i.append(y)
        all_f_i.append(f(X))
    
    optimize.minimize(
      f, x0, jac=grad, tol=tol,
      method="BFGS", callback=store 
    )
    
    return all_x_i, all_y_i, all_f_i
```
:::
::::


## Trajectory

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.7)
opt = bfgs((1.6, 1.1), f, grad)
plot_2d_traj((-1,2), (-1,2), f, traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.05)
opt = bfgs((1.6, 1.1), f, grad)
plot_2d_traj((-1,2), (-1,2), f, traj=opt)
```
:::
::::


## Rosenbrock's function 

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = bfgs((1.6, 1.1), f, grad)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = bfgs((-0.5, 0), f, grad)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::
::::


## Method: Nelder-Mead

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

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| code-line-numbers: 1,12-15
def nelder_mead(x0, f, grad, hess=None, tol=1e-8):
    all_x_i = [x0[0]]
    all_y_i = [x0[1]]
    all_f_i = [f(x0)]
    
    def store(X):
        x, y = X
        all_x_i.append(x)
        all_y_i.append(y)
        all_f_i.append(f(X))
    
    optimize.minimize(
      f, x0, tol=tol,
      method="Nelder-Mead", callback=store 
    )
    
    return all_x_i, all_y_i, all_f_i
```
:::
::::


## Nelder-Mead

::: {.center .large}
[Live Demo](http://nelder-mead.s3-website.us-east-2.amazonaws.com/)
:::

::: {.aside}
From [github.com/greg-rychlewski/nelder-mead](https://github.com/greg-rychlewski/nelder-mead)
:::


## Trajectory

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.7)
opt = nelder_mead((1.6, 1.1), f, grad)
plot_2d_traj((-1,2), (-1,2), f, traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_quad(0.05)
opt = nelder_mead((1.6, 1.1), f, grad)
plot_2d_traj((-1,2), (-1,2), f, traj=opt)
```
:::
::::


## Rosenbrock's function 

:::: {.columns .small}
::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = nelder_mead((1.6, 1.1), f, grad)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 95%
f, grad, hess = mk_rosenbrock()
opt = nelder_mead((-0.5, 0), f, grad)
plot_2d_traj((-2,2), (-2,2), f, traj=opt)
```
:::
::::


