# 2022-04-08 Optimization

## Contents

# 2022-04-08 Optimization¶

## Last time¶

Project feedback

Computing derivatives in maintainable code

Forward and reverse

Algorithmic (automatic) differentiation

Projects

## Today¶

Reflect on differentiation

Second order (Newton type) optimization

Project discussion

```
using LinearAlgebra
using Plots
default(linewidth=4, legendfontsize=12)
function vander(x, k=nothing)
if isnothing(k)
k = length(x)
end
m = length(x)
V = ones(m, k)
for j in 2:k
V[:, j] = V[:, j-1] .* x
end
V
end
function vander_chebyshev(x, n=nothing)
if isnothing(n)
n = length(x) # Square by default
end
m = length(x)
T = ones(m, n)
if n > 1
T[:, 2] = x
end
for k in 3:n
#T[:, k] = x .* T[:, k-1]
T[:, k] = 2 * x .* T[:,k-1] - T[:, k-2]
end
T
end
runge(x) = 1 / (1 + 10*x^2)
runge_noisy(x, sigma) = runge.(x) + randn(size(x)) * sigma
function grad_descent(loss, grad, c0; gamma=1e-3, tol=1e-5)
"""Minimize loss(c) via gradient descent with initial guess c0
using learning rate gamma. Declares convergence when gradient
is less than tol or after 500 steps.
"""
c = copy(c0)
chist = [copy(c)]
lhist = [loss(c)]
for it in 1:500
g = grad(c)
c -= gamma * g
push!(chist, copy(c))
push!(lhist, loss(c))
if norm(g) < tol
break
end
end
(c, hcat(chist...), lhist)
end
function diff_wp(f, x; eps=1e-8)
"""Diff using Walker and Pernice (1998) choice of step"""
h = eps * (1 + abs(x))
(f(x+h) - f(x)) / h
end
```

```
diff_wp (generic function with 1 method)
```

# Hand-coding derivatives¶

```
function f(x)
y = x
for _ in 1:2
a = y^pi
b = cos(a)
c = log(y)
y = b * c
end
y
end
f(1.9), diff_wp(f, 1.9)
```

```
(-1.5346823414986814, -34.032439961925064)
```

```
function df(x, dx)
y = x
dy = dx
for _ in 1:2
a = y ^ pi
da = pi * y^(pi - 1) * dy
b = cos(a)
db = -sin(a) * da
c = log(y)
dc = 1/y * dy
y = b * c
dy = db * c + b * dc
end
y, dy
end
df(1.9, 1)
```

```
(-1.5346823414986814, -34.03241959914048)
```

# We can go the other way¶

We can differentiate a composition \(h(g(f(x)))\) as

What we’ve done above is called “forward mode”, and amounts to placing the parentheses in the chain rule like

The expression means the same thing if we rearrange the parentheses,

which we can compute with in reverse order via

# A reverse mode example¶

```
function g(x)
a = x^pi
b = cos(a)
c = log(x)
y = b * c
y
end
(g(1.4), diff_wp(g, 1.4))
```

```
(-0.32484122107701546, -1.2559760384500684)
```

```
function gback(x, y_)
a = x^pi
b = cos(a)
c = log(x)
y = b * c
# backward pass
c_ = y_ * b
b_ = c * y_
a_ = -sin(a) * b_
x_ = 1/x * c_ + pi * x^(pi-1) * a_
x_
end
gback(1.4, 1)
```

```
-1.2559761698835525
```

# How does Zygote work?¶

```
square(x) = x^2
@code_llvm square(1.5)
```

```
; @ In[8]:1 within `square`
define double @julia_square_5092(double %0) #0 {
top:
; ┌ @ intfuncs.jl:312 within `literal_pow`
; │┌ @ float.jl:405 within `*`
%1 = fmul double %0, %0
; └└
ret double %1
}
```

```
@code_llvm Zygote.gradient(square, 1.5)
```

```
; @ /home/jed/.julia/packages/Zygote/H6vD3/src/compiler/interface.jl:74 within `gradient`
define [1 x double] @julia_gradient_5209(double %0) #0 {
top:
; @ /home/jed/.julia/packages/Zygote/H6vD3/src/compiler/interface.jl:76 within `gradient`
; ┌ @ /home/jed/.julia/packages/Zygote/H6vD3/src/compiler/interface.jl:41 within `#56`
; │┌ @ In[8]:1 within `Pullback`
; ││┌ @ /home/jed/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67 within `#1822#back`
; │││┌ @ /home/jed/.julia/packages/Zygote/H6vD3/src/lib/number.jl:6 within `#254`
; ││││┌ @ promotion.jl:380 within `*` @ float.jl:405
%1 = fmul double %0, 2.000000e+00
; └└└└└
; @ /home/jed/.julia/packages/Zygote/H6vD3/src/compiler/interface.jl:77 within `gradient`
%.fca.0.insert = insertvalue [1 x double] zeroinitializer, double %1, 0
ret [1 x double] %.fca.0.insert
}
```

# Kinds of algorithmic differentation¶

Source transformation: Fortran code in, Fortran code out

Duplicates compiler features, usually incomplete language coverage

Produces efficient code

Operator overloading: C++ types

Hard to vectorize

Loops are effectively unrolled/inefficient

Just-in-time compilation: tightly coupled with compiler

JIT lag

Needs dynamic language features (JAX) or tight integration with compiler (Zygote, Enzyme)

Some sharp bits

# Forward or reverse?¶

It all depends on the shape.

One input, many outputs: use forward mode

“One input” can be looking in one direction

Many inputs, one output: use reverse mode

Will need to traverse execution backwards (“tape”)

Hierarchical checkpointing

Small intermediate results?

Break into pieces and use forward/reverse for the smaller pieces.

About square? Forward is usually a bit more efficient.

# Can you differentiate an algorithm?¶

Input \(c\), output \(x\) such that \(f(x,c) = 0\)

Input \(A\), output \(\lambda\) such that \(A x = \lambda x\) for some nonzero vector \(x\)

Input

`buffer`

, output`sha256(buffer)`

# Ill-conditioned optimization¶

Gradient of \(L\) requires the Jacobian \(J\) of the model \(f\).

We can solve \(g(c) = 0\) using a Newton method

The Hessian requires the second derivative of \(f\), which can cause problems.

# Newton-like methods for optimization¶

Solve

Update \(c \gets c + \gamma \delta c\) using a line search or trust region.

Gauss-Newton: \(H = J^T J\)

Levenberg-Marquardt: \(H = J^T J + \alpha I\)

# Outlook¶

The optimization problem can be solved using a Newton method. It can be onerous to implement the needed derivatives.

The Gauss-Newton method is often more practical than Newton while being faster than gradient descent, though it lacks robustness.

The Levenberg-Marquardt method provides a sort of middle-ground between Gauss-Newton and gradient descent.

Many globalization techniques are used for models that possess many local minima.

One pervasive approach is stochastic gradient descent, where small batches (e.g., 1 or 10 or 20) are selected randomly from the corpus of observations (500 in our current example), and a step of gradient descent is applied to that reduced set of observations. This helps to escape saddle points and weak local minima.

Among expressive models \(f(x,c)\), some may converge much more easily than others. Having a good optimization algorithm is essential for nonlinear regression with complicated models, especially those with many parameters \(c\).

Classification is a very similar problem to regression, but the observations \(y\) are discrete, thus

models \(f(x,c)\) must have discrete output

the least squares loss function is not appropriate.

# Team project workshopping¶

Identify one or more software packages and execution environments that will be used.

Identify one or more stakeholders in this problem space.

Identify one or more decisions that depend on the use or behavior of the methods or software.

What will you do?

Identify a test you can conduct to inform the decision. Think about the stakeholders and context.

Create something new using the software and reflect on how it works.

Make a contribution to improve the software or community (diagnostics, decumentation, efficiency). (Please share your ideas with me before choosing this.)

Create a thread in

`#projects`

summarizing the above and naming your group members.