---
title: Optimization and Its Algorithms
author: 36-465/665, Spring 2021
date: 11 March 2021 (Lecture 11)
output: slidy_presentation
bibliography: locusts.bib
---

\[
\newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right)}
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]}
\newcommand{\Risk}{r}
\newcommand{\EmpRisk}{\hat{r}}
\newcommand{\Loss}{\ell}
\newcommand{\OptimalStrategy}{\sigma}
\DeclareMathOperator*{\argmin}{argmin}
\newcommand{\ModelClass}{S}
\newcommand{\OptimalModel}{s^*}
\DeclareMathOperator{\tr}{tr}
\newcommand{\Indicator}[1]{\mathbb{1}\left\{ #1 \right\}}
\newcommand{\myexp}[1]{\exp{\left( #1 \right)}}
\newcommand{\eqdist}{\stackrel{d}{=}}
\newcommand{\Rademacher}{\mathcal{R}}
\newcommand{\EmpRademacher}{\hat{\Rademacher}}
\newcommand{\Growth}{\Pi}
\newcommand{\VCD}{\mathrm{VCdim}}
\newcommand{\OptDomain}{\Theta}
\newcommand{\OptDim}{p}
\newcommand{\optimand}{\theta}
\newcommand{\altoptimand}{\optimand^{\prime}}
\newcommand{\ObjFunc}{M}
\newcommand{\outputoptimand}{\optimand_{\mathrm{out}}}
\newcommand{\Hessian}{\mathbf{h}}
\]


## Previously

- Empirical risk of strategy $s$ is $\EmpRisk(s) = n^{-1}\sum_{i=1}^{n}{\Loss(Y_i, s(X_i))}$
- We want to find the empirical risk minimizer $\hat{s}$,
\[
\hat{s} \equiv \argmin_{s \in \ModelClass}{\EmpRisk(s)}
\]
- We've used properties of the model class $\ModelClass$ to control generalization error
- We've used how $\ModelClass$ behaves under _this_ distribution of data to control generalization error (Rademacher complexity)
- What about using properties of the actual fitting process?
- We're now going to open up the black box of $\argmin$


## Optimization basics

- The function we're trying to optimize is the **objective function**, let's say $\ObjFunc$ today
- The argument to $\ObjFunc$ is (say) $\optimand$
    + Some people call this the **optimand**
- The possible values of $\optimand$ is $\OptDomain$, the **domain** or **feasible set**, whose dimension is (say) $\OptDim$
- **Optimization** can be minimization or maximization, as we like; we'll stick with minimizing


## Local vs. global minima

- $\optimand$ is a **global minimum** when $\altoptimand \neq \optimand$ $\Rightarrow$ $\ObjFunc(\altoptimand) \geq \ObjFunc(\optimand)$
    + Not necessarily unique!
- $\optimand$ is a **local minimum** when $\ObjFunc(\optimand) \leq \ObjFunc(\altoptimand)$ whenever $\altoptimand$ is close enough to $\optimand$
    + Every global minimum is also a local minimum
    + If there's only one local minimum anywhere, it's the global minimum
- Lots of local minima tend to make it harder to find the global minimum
	

## "The" minimum: value vs. location

- If $\optimand^*$ is a global minimum, then $\ObjFunc(\optimand^*)$ is the **value** of the minimum or minimal value, in symbols
\[
\min_{\optimand \in \OptDomain}{\ObjFunc(\optimand)}
\]
- But $\optimand^*$ itself is the **location** of the global minimum, in symbols
\[
\argmin_{\optimand \in \OptDomain}{\ObjFunc(\optimand)}
\]
- Example: the minimal value of $(x-1)^2$ is 0, but the location of the minimum is $x=1$
- Both value and location can change with $\OptDomain$; this will be important when we look at constraints later


## Finding the minimum: optimization algorithms

- An optimization algorithm starts from $\ObjFunc$ and $\OptDomain$, and possibly an initial guess $\optimand^{(0)}$, and returns an **approximation** to $\argmin_{\optimand \in \OptDomain}{\ObjFunc(\optimand)}$,  call it $\outputoptimand$
    + Really the algorithm takes in **representations** of all these things, and sometimes there are interesting issues with better or worse representations, amount of detail in the representation [@Traub-Werschulz-complexity-and-info], etc.
- We usually measure the approximation by difference in the _value_, not the location: the algorithm gets $\epsilon$-close when
\[
\ObjFunc(\outputoptimand) \leq \epsilon + \min_{\optimand \in \OptDomain}{\ObjFunc(\optimand)}
\]
- Often, the longer we let the algorithm run, the better the approximation
    + How many steps does the algorithm need to get $\epsilon$-close to the optium?  $O(1/\epsilon)$ or $O(\epsilon^{-d})$ would be polynomial, $O(\log{1/\epsilon})$ would be logarithmic, $\myexp{O(1/\epsilon)}$ would be exponential (and bad)


## So how do we build an optimization algorithm anyway?

- Start with calculus, assume $\ObjFunc$ has as many derivatives as we need
    + All derivatives are with respect to $\optimand$
- At an (interior!) local minimum, or local maximum, or inflection point,
\[
\nabla \ObjFunc(\optimand) = 0
\]
    + $\ObjFunc$ is flat at the minimum
    + In 1D, this is just $d\ObjFunc/d\optimand = 0$
    + This is called the **first-order condition** (or conditions)
- At an (interior) local _minimum_, but not any other sort of extremal point, for any $\vec{v} \neq 0$,
\[
\vec{v} \cdot \nabla\nabla \ObjFunc(\optimand) \vec{v} \geq 0
\]
    + and it's positive in at least some directions!
    + Remember $\nabla\nabla \ObjFunc$ is the matrix of 2nd partial derivatives, or **Hessian**, so I will also write $\Hessian$
    + Jargon: this condition is what's meant when we say the Hessian is **non-negative-definite**, in symbols $\Hessian(\theta) \succeq 0$
	    + **Positive-definite**, $\Hessian(\theta) \succ 0$, means $\vec{v} c\dot \Hessian \vec{v} > 0$ for all $\vec{v} \neq 0$; this implies $\theta$ is an _isolated_ local minimum
    + In words: moving away from the local minimum in any direction can only increase the objective function
    + In 1D this is just $\frac{d^2\ObjFunc}{d\optimand^2} > 0$
- $\nabla \ObjFunc(\optimand)=0$ is the **first-order condition** and $\Hessian(\optimand)  \succeq 0$ is the **second-order condition**


## Optimizing by equation-solving

- One approach: use the first-order condition to get a system of equations
\[
\nabla \ObjFunc(\optimand) = 0
\]
    + One equation per coordinate of $\optimand$
    + When $\ObjFunc$ is empirical risk $\EmpRisk$, sometimes called the **estimting equations** or even **normal equations**
- Solve the system of equations for $\optimand$
- If there's more than one solution, check the second-order conditions
- This is what you did for estimating linear models by ordinary least squares, or even weighted least squares


## Pros and cons of the solve-the-equations approach

> - You need to _set up_ the system of equations, and often finding $\nabla \ObjFunc$ would itself be a pain
>    + Numerical differentiation is a thing, however
> - You need to solve a system of equations: good if there are good **solvers** for that type of system of equations, not so good otherwise
>    + 200+ years of work have given us _very_ good solvers for linear systems
>        * For linear systems, even very old-fashioned methods that go back to Gauss around 1800 get $\epsilon$ approximations with $O(\log{1/\epsilon})$ iterations
>    + General-purpose nonlinear equation-solving is still much harder
>        * sometimes works by using Taylor expansion to linearize
>        * sometimes works by turning the solve-the-equations into "minimize the difference between the left and the right hand side of the equation"


## Go back to the calculus

- Start with a guess $\optimand^{(0)}$
- Find $\nabla \ObjFunc(\optimand^{(0)})$
- Move in the _opposite_ direction:
\[
\optimand^{(1)} = \optimand^{(0)} - a_0 \nabla \ObjFunc(\optimand^{(0)})
\]
- Repeat:
\[
\optimand^{(t+1)} = \optimand^{(t)} - a_t \nabla \ObjFunc(\optimand^{(t)})
\]
- Notice: a local optimum will be a fixed point!
- Issue: how big are the **step sizes** $a_t$?
    + Confusingly also called the **learning rate**
	
## Constant-step-size gradient descent

- Inputs: objective function $\ObjFunc$, step size $a$, initial guess $\optimand^{(0)}$
- While ((not too tired) and (making adequate progress))
    + Find $\nabla \ObjFunc(\optimand^{(t)})$
    + Set $\optimand^{(t+1)} \leftarrow \optimand^{(t)} - a \nabla \ObjFunc(\optimand^{(t)})$
- Return final $\optimand$

- "not too tired": Set a maximum number of iterations
- "making adequate progress": $ObjFunc$ isn't changing by too little
to bother with, $\optimand$ isn't changing by too little to bother with,
$\nabla \ObjFunc$ isn't too close to zero


## Constant-step-size gradient descent

- Pick an $a > 0$ that's small and use it at each step
- Each iteration of gradient descent takes $O(\OptDim))$ operations
    + Find $\OptDim$ derivatives, multiply by $a$, add to $\optimand^{(t-1)}$
- If $\ObjFunc$ is nice, $\optimand^{(t)}$ is an $\epsilon$-approximation of the optimum after $t=O(\epsilon^{-2})$ iterations
    + i.e. at that point $\ObjFunc(\optimand^{(t)}) \leq \epsilon+\min{\ObjFunc(\optimand)}$
	+ "Nice" here means: convex and second-differentiable
- If $\ObjFunc$ is _very_ nice, $\optimand^{(t)}$ is an $\epsilon$-approximation after only $t=O(\log{1/\epsilon})$ iterations
    + "nice" plus _strictly_ convex
	

## Gradient descent is basic, but powerful

- Gradient descent works well when there's a single global minimum, no flat parts to the function, and the step size is small enough to not over-shoot or zig-zag
- It's actually been re-invented a number of times under different names
     + e.g., "back-propagation" [@Rumehlart-Hinton-Williams-introduce-backprop]
 - It's the work-horse for large-scale industrial applications in modern machine learning
    + especially as **stochastic gradient descent**
- It's still a bit mysterious why it works so well for those applications, which actually have lots of local minima!


## Estimation error vs. optimization error

- Remember our approximation error vs. estimation error decomposition:
\begin{eqnarray}
\Risk(\hat{s}) & = & \Risk(\OptimalStrategy) + (\Risk(\OptimalModel) - \Risk(\OptimalStrategy)) + (\Risk(\hat{s}) - \Risk(\OptimalModel))\\
& = & (\text{true minimum risk}) + (\text{approximation error from limited strategy set})\\
& & + (\text{estimation error from not knowing the best-in-class set})
\end{eqnarray}
- Now we don't even have $\hat{s} = \argmin{\EmpRisk(s)}$, we have $\hat{s}_{\mathrm{out}}$, the output of some algorithm
\begin{eqnarray}
\Risk(\hat{s}_{\mathrm{out}}) & =  & \Risk(\OptimalStrategy) + (\Risk(\OptimalModel) - \Risk(\OptimalStrategy)) + (\Risk(\hat{s}) - \Risk(\OptimalModel)) + (\Risk(\hat{s}_{\mathrm{out}}) - \Risk(\hat{s}))\\
&= & (\text{optimal risk}) + (\text{approximation error}) + (\text{estimation error}) + (\text{optimization error})
\end{eqnarray}
- Optimization error $\approx$ what I've been calling $\epsilon$
    + only $\approx$ because of $\Risk$ vs. $\EmpRisk$ issue


## Estimation error vs. optimization error (2)

\[
\text{risk} = \text{minimal risk} + \text{approximation error} + \text{estimation error} + \text{optimization error}
\]

- Minimal risk and approximation error don't change with $n$ or with how we optimize
- Estimation error shrinks with $n$: for large $n$, typically $O(\OptDim/n)$
    + Possibly more slowly converging in $n$ for some families, or if $\OptDim$ grows with $n$
    + VC theory tells us estimation error can be $O(\frac{\log{n/d}}{n/d})$
- Optimization error shrinks as we do more computational work
- There's no point to making the optimization error much smaller than the estimation error
    + More exactly: lots of work for little real benefit
- So: don't try to make the optimization error much smaller than $O(\OptDim/n)$

> Don't bother optimizing more precisely than the noise in the data will support


## Beyond gradient descent: Newton's method

- Needing to pick the step-size $a_t$ is annoying
- We'd like to take big steps, but $\nabla \ObjFunc$ is a local quantity and might be mis-leading far away
- $\Rightarrow$ We'd like to take bigger steps when the gradient doesn't change much
- This is **Newton's method**:
\[
\optimand^{(t+1)} = \optimand^{(t)} - \left(\Hessian(\optimand^{(t)}) \right)^{-1} \nabla \ObjFunc(\optimand^{(t)})
\]
    + One route to this: pretend $\ObjFunc$ is quadratic, as justified by a Taylor expansion around the true minimum
- This is like gradient descent, but using the inverse Hessian to give the step size
    + And possibly a bit of rotation away from the gradient
	
## Pros of Newton's method

- Adaptively-chosen step size means harder to get zig-zags, over-shooting, etc.
- Need $O(\epsilon^{-2})$ steps to get an $\epsilon$ approximation to the minimum for nice functions
- For very nice functions, only need $O(\log{(\log{(1/\epsilon)})})$ iterations
- Generally needs many fewer iterations than gradient descent


## Cons of Newton's method

- Hopeless if the Hessian doesn't exist or isn't invertible
- Need to take $O(p^2)$ second derivatives _and_ $p$ first derivatives, total $O(p^2)$
- Need to find $\optimand^{(t+1)}$
    + Seems straightforward, it's $\optimand^{(t+1)} = \optimand^{(t)} - \left(\Hessian(\optimand^{(t)}) \right)^{-1} \nabla \ObjFunc(\optimand^{(t)})$
    + But inverting a $[p\times p]$ matrix takes $O(p^3)$ operations in general, so this would be an $O(p^3)$ step
- Alternative: solve $\Hessian \optimand^{(t+1)} = \Hessian \optimand^{(t)} - \nabla \ObjFunc(\optimand^{(t)})$ for $\optimand^{(t+1)}$ for the unknown $\optimand^{(t+1)}$
    + (Take the basic update equation for Newton's method and multiply both sides by $\Hessian$ from the left)
    + Solving a system of $p$ linear equations can be done in $O(p^2)$ time for a _particular_ right-hand side
    + Lots of variants to use approximate Hessians rather than the full deal (BFGS, built in to R's `optim()`, is one of these)
- So each iteration is $O(p^2)$, much slower than gradient descent's $O(p)$
    + $O(p^2)$ to get Hessian and gradient plus $O(p^2)$ to solve for update $=O(p^2)$


## Gradient methods with big data

\[
\EmpRisk(\theta) = \frac{1}{n}\sum_{i=1}^{n}{\Loss(y_i, s(x_i; \theta))}
\]

- Getting a value of $\EmpRisk$ at a particular $\theta$ is $O(n)$, getting $\nabla \EmpRisk$ is $O(np)$, getting $\Hessian$ is $O(np^2)$
    + And that's assuming calculating $s(x_i; \theta)$ doesn't slow down with $n$
- Maybe OK when $n=100$ or $n=10^4$, but with $n=10^9$ or $n=10^{12}$, we really don't know which way to move


## Gradient methods with big data (2)

![](sarcastic-gradient-descent.png)

[(src)](https://twitter.com/ml_hipster/status/226448024224796673)



## A way out: sampling is an unbiased estimate

- Pick _one_ data point $I$ at random, uniform on $1:n$
- $\Loss(y_I, s(x_I; \theta))$ is random, but
\[
\Expect{\Loss(y_I, s(x_I; \theta))} = \EmpRisk(\theta)
\]
- Re-brand $\Loss(y_I, s(x_I; \theta))$ as $\EmpRisk_{I}(\theta)$
\begin{eqnarray}
\Expect{\EmpRisk_{I}(\theta)} & = & \EmpRisk(\theta)\\
\Expect{\nabla \EmpRisk_{I}(\theta)} & = & \nabla \EmpRisk(\theta)\\
\Expect{\nabla \nabla \EmpRisk_{I}(\theta)} & = & \Hessian(\theta)
\end{eqnarray}
- $\Rightarrow$ Don't optimize with all the data, optimize with random samples

## Stochastic gradient descent

- Draw _lots_ of random one-point samples and let their noise cancel out:

0. Start with initial guess $\optimand^{(0)}$, adjustment rate $a$
1. While (not too tired) and (making adequate progress))
    a. At $t^{\mathrm{th}}$ iteration, pick random $I$ uniformly on $1:n$
	b. Set $\optimand^{(t+1)} \leftarrow \optimand^{(t)} - \frac{a}{t}\nabla \EmpRisk_{I}(\optimand^{(t)})$
2. Return final $\optimand$

- Shrinking step-sizes by $1/t$ ensures noise in each gradient dies down

## Stochastic gradient descent (2)

- Tons of variants:
    + Put the data points $1:n$ in a random order and then cycle through them
    + Don't check the "making adequate progress" condition too often
    + Adjust the $1/t$ step-size to some other function
    + **Stochastic Newton's method**: Use the sample to also calculate the Hessian and take a Newton's method step
    + **Mini-batch**: Sample a _few_ of random data points at once
    + Mini-batch stochastic Newton's method, etc.
	
	
## Pros and cons of stochastic gradient methods

- Pro: Each iteration is (or at least constant in $n$)
- Pro: Never need to hold all the data in memory at once
- Pro: Does converge eventually (at least if the non-stochastic method would)
- Cons: sampling noise increases optimization error
    + That is: more iterations to come within the same $\epsilon$ of the optimum as non-stochastic GD or Newton
- Over-all pro: often low computational cost to make the optimzation error small compared to the estimation error


## More optimization algorithms

- Ones which play more tricks with derivatives than just gradient descent and Newton ("conjugate gradient", etc., etc.)
- Ones which avoid derivatives ("simplex" or "Nelder-Mead")
- Ones which avoid derivatives and try random changes ("simulated annealing")
- Ones which use natural-selection-with-random-variation to evolve a whole population of approximate optima ("genetic algorithms")


## Why are there so many different optimization algorithms?

> - "Come up with a new algorithm" is a way to make a mark...
> - No one algorithm works well on _every_ problem
>     + Sometimes obvious: don't use Newton's method if $\OptDomain$ is discrete
> - Fundamental limit: no algorithm is _universally_ better than others on _every_ problem (**no free lunch theorem** of @Wolpert-Macready-no-free-lunch)
>    + For every problem where your favorite algorithm does better than mine, I can design a new problem where my algorithm leads by just as much [@Culberson-no-free-lunch]
> - We need to know _something_ about the problem to select a good optimizer


## Summing up

> - With real data and real computers, finding the empirical-risk-minimizer means using an algorithm to solve an optimization problem
> - These algorithms almost never give the _exact_ optimum but just an approximation
> - Usually, the longer an algorithm is allowed to work, the closer it can get to the true optimum
> - This adds optimization error on to estimation error
> - For many statistical learning problems, gradient descent and Newton's method work really well
>     + With sampling to make them more computationally efficient for big data
> - Don't bother reducing the optimization error much beyond the estimation error
> - No algorithm is best for all problems


## References




