# Automatic differentiation and nonliner optimization

This notebook is based on:

Arthur Delarue's notebook from COS 2018: https://github.com/PhilChodrow/cos_2018/blob/master/7_nonlinear_and_integer/Automatic%20Differentiation%20and%20Dual%20Numbers%20-%20Solutions.ipynb
material from https://nbviewer.jupyter.org/github/JuliaOpt/JuMPTutorials.jl/blob/master/notebook/using_JuMP/nonlinear_modelling.ipynb contributed to by Arpit Bhatia.

## Computing derivatives for nonlinear optimization using automatic differentiation
Consider a general constrained nonlinear optimization problem:

$$
\begin{align}
\min \quad & f(x) \\
\text{s.t.} \quad & g(x) = 0 \\
& h(x) \leq 0.
\end{align}
$$

where $f : \mathbb{R}^n \to \mathbb{R}, g : \mathbb{R}^n \to \mathbb{R}^r$, and $h: \mathbb{R}^n \to \mathbb{R}^s$.

When $f$ and $h$ are convex and $g$ is affine, we can hope for a globally optimal solution, otherwise typically we can only ask for a locally optimal solution.

What approaches can we use to solve this?

- When $r=0$ and $s = 0$ (unconstrained), and $f$ differentiable, most classical approach is [gradient descent](http://en.wikipedia.org/wiki/Gradient_descent), also fancier methods like [Newton's method](http://en.wikipedia.org/wiki/Newton%27s_method) and quasi-newton methods like [BFGS](http://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm).
- When $f$ differentiable and $g$ and $h$ linear, [gradient projection](http://neos-guide.org/content/gradient-projection-methods),
- When $f$, $g$, and $h$ twice differentiable, [interior-point methods](http://en.wikipedia.org/wiki/Interior_point_method), [sequential quadratic programming](http://www.neos-guide.org/content/sequential-quadratic-programming)
- When derivatives not available, [derivative-free optimization](http://rd.springer.com/article/10.1007/s10898-012-9951-y)

This is not meant to be an exhaustive list, see http://plato.asu.edu/sub/nlores.html#general and http://www.neos-guide.org/content/nonlinear-programming for more details.

## Dual Numbers


Consider numbers of the form $x + y\epsilon$ with $x,y \in \mathbb{R}$. We *define* $\epsilon^2 = 0$, so,
$$
(x_1 + y_1\epsilon)(x_2+y_2\epsilon) = x_1x_2 + (x_1y_2 + x_2y_1)\epsilon.
$$

These are called the *dual numbers*. Think of $\epsilon$ as an infinitesimal perturbation (you've probably seen hand-wavy algebra using $(dx)^2 = 0$ when computing integrals - this is the same idea).

How can these funky *dual numbers* help us?
If we are given an infinitely differentiable function in Taylor expanded form
$$
f(x) = \sum_{k=0}^{\infty} \frac{f^{(k)}(a)}{k!} (x-a)^k
$$
it follows that 
$$
f(x+y\epsilon) = \sum_{k=0}^{\infty} \frac{f^{(k)}(a)}{k!} (x-a+y\epsilon)^k = \sum_{k=0}^{\infty} \frac{f^{(k)}(a)}{k!} (x-a)^k + y\epsilon\sum_{k=0}^{\infty} \frac{f^{(k)}(a)}{k!}\binom{k}{1} (x-a)^{k-1} = f(x) + yf'(x)\epsilon,
$$

Let's unpack what's going on here. We started with a function $f : \mathbb{R} \to \mathbb{R}$. Dual numbers are *not* real numbers, so it doesn't even make sense to ask for the value $f(x+y\epsilon)$ given $x+y\epsilon \in \mathbb{D}$ (the set of dual numbers). But anyway we plugged the dual number into the Taylor expansion, and by using the algebra rule $\epsilon^2 = 0$ we found that $f(x+y\epsilon)$ must be equal to $f(x) + yf'(x)\epsilon$ if we use the Taylor expansion as the definition of $f : \mathbb{D} \to \mathbb{R}$.


Alternatively, for any once differentiable function $f : \mathbb{R} \to \mathbb{R}$, we can define its extension to the dual numbers as$$
f(x+y\epsilon) = f(x) + yf'(x)\epsilon.
$$This is essentially equivalent to the previous definition.

Let's verify a very basic property, the chain rule, using this definition.

Suppose $h(x) = f(g(x))$. Then,$$
h(x+y\epsilon) = f(g(x+y\epsilon)) = f(g(x) + yg'(x)\epsilon) = f(g(x)) + y \underbrace{g'(x)f'(g(x))}_{=\text{chain rule}} \epsilon = h(x) + yh'(x)\epsilon.
$$

Maybe that's not too surprising, but it's actually a quite useful observation.

## Implementation

We will use the package `ForwardDiff`, but there are others: e.g. `ReverseDiff`, `Zygote`


```julia
using ForwardDiff
```


```julia
d = ForwardDiff.Dual(3.0, 4.0)
```




    Dual{Nothing}(3.0,4.0)




```julia
typeof(d)
```




    ForwardDiff.Dual{Nothing, Float64, 1}




```julia
typeof(d)
```




    ForwardDiff.Dual{Nothing, Float64, 1}




```julia
ForwardDiff.partials(d)
```




    1-element ForwardDiff.Partials{1, Float64}:
     4.0




```julia
@which d + ForwardDiff.Dual(3.0,4.0)
```




+(x::<b>ForwardDiff.Dual{Txy, V, N} where N where V</b>, y::<b>ForwardDiff.Dual{Txy, V, N} where N where V</b>)<i> where Txy</i> in ForwardDiff at <a href="file://C:/Users/lkape/.julia/packages/ForwardDiff/qTmqf/src/dual.jl" target="_blank">C:\Users\lkape\.julia\packages\ForwardDiff\qTmqf\src\dual.jl:134</a>



The effect should be:

`+(z::Dual, w::Dual) = dual(real(z)+real(w), epsilon(z)+epsilon(w))`


```julia
@which ForwardDiff.Dual(2.0,2.0) * ForwardDiff.Dual(3.0,4.0)
```




*(x::<b>ForwardDiff.Dual{Txy, V, N} where N where V</b>, y::<b>ForwardDiff.Dual{Txy, V, N} where N where V</b>)<i> where Txy</i> in ForwardDiff at <a href="file://C:/Users/lkape/.julia/packages/ForwardDiff/qTmqf/src/dual.jl" target="_blank">C:\Users\lkape\.julia\packages\ForwardDiff\qTmqf\src\dual.jl:134</a>



The effect should be:

`*(z::Dual, w::Dual) = dual(real(z)*real(w), epsilon(z)*real(w)+real(z)*epsilon(w))`

## How general is it?


```julia
function squareroot(x)
    z = x # Initial starting point
    while abs(z * z - x) > 1e-13
        z = z - (z * z - x) / (2z)
    end
    return z
end
```




    squareroot (generic function with 1 method)




```julia
squareroot(100)
```




    10.0



Can we differentiate this code? Yes!


```julia
d = squareroot(ForwardDiff.Dual(100.0, 1.0))
```




    Dual{Nothing}(10.0,0.049999999999999996)




```julia
1 / (2 * sqrt(100)) # The exact derivative
```




    0.05



If we were to do this ourselves:


```julia
function squareroot(x)
    z = x # Initial starting point
    z_prime = 1
    while abs(z * z - x) > 1e-13
        z = z - (z * z - x) / (2z)
        # apply the quotient rule to the above:
        z_prime = z_prime - ((2z * z_prime - 1) * 2z - (z * z - x) * 2z_prime) / (4 * z ^ 2)
    end
    return (z, z_prime)
end
squareroot(100)
```




    (10.0, 0.05)



## In optimization

The Gaussian likelihood is -

$$
L(\theta | \mathbf{x})=L\left(\theta_{1}, \ldots, \theta_{k} | x_{1}, \ldots, x_{n}\right)=\prod_{i=1}^{n} f\left(x_{i} | \theta_{1}, \ldots, \theta_{k}\right)
$$

Suppose we want to maximize log-likelihood. 
We can use `@NLobjective` to do this in JuMP and a nonlinear solver e.g. IPOPT.

IPOPT will use derivative information.
We are allowed to plug-in log-likelihood directly into `@NLobjective` because it so happens that our log-likelihood calculation also uses functions with derivatives that are already defined in `Calculus.jl`, which JuMP will use.

The following example is taken from https://nbviewer.jupyter.org/github/JuliaOpt/JuMPTutorials.jl/blob/master/notebook/using_JuMP/nonlinear_modelling.ipynb.

### Exercise

Use `@NLobjective` so that we minimize the log-likelihood.
The likelihood is given by:
$$
\left( 2 \pi \sigma^2 \right)^{-n/2} \exp \left( - \sum_{i=1}^n (x_i - \mu)^2 / (2 \sigma^2) \right)
$$
So the log-likelihood is:
$$
-\frac{n}{2} \log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_{i=1}^n (x_i - \mu)^2.
$$


```julia
using JuMP
using Ipopt

function maximize_log_likelihood(x::Vector{Float64})
    n = length(x)
    mle = Model(Ipopt.Optimizer)
    set_parameters(mle, "print_level" => 0)

    @variable(mle, μ, start = 0.5)
    @variable(mle, σ >= 0.0, start = 0.5)
    
    @NLobjective(mle, Max, log((2 * π * σ^2)^(-n / 2) * exp(-(sum((x[i] - μ)^2 for i in 1:n) / (2 * σ^2)))))
    
     optimize!(mle)
    return (value(μ), value(σ))
end
```




    maximize_log_likelihood (generic function with 1 method)




```julia
# test your function
using Random
Random.seed!(1)
println(maximize_log_likelihood(randn(10)))
println(maximize_log_likelihood(randn(50)))
println(maximize_log_likelihood(randn(100)))
```

    
    ******************************************************************************
    This program contains Ipopt, a library for large-scale nonlinear optimization.
     Ipopt is released as open source code under the Eclipse Public License (EPL).
             For more information visit http://projects.coin-or.org/Ipopt
    ******************************************************************************
    
    (0.053307729269647004, 1.1062693505952153)
    (-0.0932713652401274, 0.9799742023885423)
    (-0.18294427514125225, 1.0528182497375815)


What if we had a more involved function e.g. a function with a loop, output from a simulation, output from an ODE, or anything that is too tedious to work out derivatives for manually?

We can register user-defined functions. JuMP will use AD to get derivative information.


```julia
Random.seed!(1)
x = randn(10);
```


```julia
function loglikelihood(μ::Float64)
    # assume we have σ = 1
    n = length(x)
    # make sure to work in logspace so that numbers don't "explode"
    return -n / 2 * log(2 * π) - sum((x[i] - μ)^2 for i in 1:n) / 2
end
```




    loglikelihood (generic function with 1 method)




```julia
grad = ForwardDiff.derivative(loglikelihood, 0.78) # doesn't work
```


```julia
# we typed our function too strictly, let's try again
function loglikelihood(μ)
    # assume we have σ = 1
    n = length(x)
    return -n / 2 * log(2 * π) - sum((x[i] - μ)^2 for i in 1:n) / 2
end
```




    loglikelihood (generic function with 2 methods)




```julia
grad = ForwardDiff.derivative(loglikelihood, 0.78)
```




    -7.26692270730353




```julia
mle = Model(Ipopt.Optimizer)
set_parameters(mle, "print_level" => 0)
register(mle, :loglikelihood, 1, loglikelihood, autodiff = true)
@variable(mle, μ, start = 1.0)
@NLobjective(mle, Max, loglikelihood(μ))
optimize!(mle)
println(termination_status(mle))
value(μ)
# previously mean was 0.053307729269647004
```

    LOCALLY_SOLVED





    0.05330772926964722




```julia

```

78584d7f89a4d694c21f26c2d2e3eaef62f7f6df
ipynb
Jupyter Notebook
6+7_julia_and_jump/.ipynb_checkpoints/AD and nonlinear complete-checkpoint.ipynb
adelarue/cos_2021
672f32637bd528b700e0e0ba099c17d4b3c26341
MIT
2020-11-14T16:12:49.000Z
2022-01-03T16:57:02.000Z
6+7_julia_and_jump/.ipynb_checkpoints/AD and nonlinear complete-checkpoint.ipynb
adelarue/cos_2021
672f32637bd528b700e0e0ba099c17d4b3c26341
MIT
6+7_julia_and_jump/.ipynb_checkpoints/AD and nonlinear complete-checkpoint.ipynb
adelarue/cos_2021
672f32637bd528b700e0e0ba099c17d4b3c26341
MIT
2021-01-28T01:35:27.000Z
2021-01-28T01:35:27.000Z
Qwen/Qwen-72B
1. YES
2. YES