# Optim user-supplied gradient and hessian MLE

**URL:** <https://discourse.julialang.org/t/optim-user-supplied-gradient-and-hessian-mle/76848>\
**Category:** Optimization (Mathematical)\
**Tags:** optim\
**Created:** [February 21, 2022, 3:14pm UTC](https://discourse.julialang.org/t/optim-user-supplied-gradient-and-hessian-mle/76848 "2022-02-21T15:14:54Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![luischanci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luischanci/32/29503_2.png) [@luischanci](https://discourse.julialang.org/u/luischanci)\
**Post date:** [February 21, 2022, 3:14pm UTC](https://discourse.julialang.org/t/optim-user-supplied-gradient-and-hessian-mle/76848/1 "2022-02-21T15:14:54Z")

</div>

I’m moving my MLE class notes to incorporate Julia, so I’m exploring how to incorporate the analytical gradient and hessian into an optimization problem using Optim. I created a toy model to check the commands first.  
However, although the code works well using only the objective function and ForwardDiff, it fails when using the gradient and/or hessian.  
Does anyone know what might be going on?  
Thanks!

Note: I already read about the potential (limited) gains obtained from incorporating the hessian, but still wanna try it.

```julia
using Optim, ForwardDiff
using LinearAlgebra, Random
Random.seed!(12); # Fix random seed generator for reproducibility

n = 500 # Number of observations
nvar = 2 # Number of variables
β = ones(nvar) * 3.71 # True coefficients
x = [ones(n) randn(n, nvar - 1)] # X matrix of explanatory variables plus constant
ε = randn(n) * 2 # Error variance
y = x * β + ε; # Generate Data

function Loglik(N,Y,X,Γ)
    σ² = Γ[end,1]
    ep = (Y - X*Γ[1:(length(Γ)-1)])
    LL = -(N/2)*log(2π) - (N/2)*log(σ²) - (ep'*ep)/(2*σ²)
    return -LL
end
function ▽f(N,Y,X,▽,Γ)
    σ² = Γ[end,1]
    k = (length(Γ)-1)
    ep = (Y - X*Γ[1:k])
    T = promote_type(eltype(X),eltype(Γ))
    ▽ .= T(0)
    ▽[1:k] = (1/σ²)*(X'*ep)
    ▽[end] = - N/(2*σ²) + (ep'*ep)/(2*(√σ²)^4)
    return ▽
end
function ▽²f(N,Y,X,H,Γ)
    σ² = Γ[end,1]
    k = (length(Γ)-1)
    ep = (Y - X*Γ[1:k])
    H[1:k,1:k] = -(1/σ²)*(X'X);
    H[1:k,end] = -(1/(σ²)^2)*(X'ep);
    H[end,1:k] = -((1/(σ²)^2)*(X'ep))';
    H[end,end] = (N/(2*(σ²)^2)) - (1/(σ²)^3)*(ep'*ep)
    return H
end

f(α) = Loglik(n,y,x,α)
g!(G,α) = ▽f(n,y,x,G,α)
h!(H,α) = ▽²f(n,y,x,H,α)
Guess = ones(size(x)[2]+1);

# prelim check that grad and hess are ok:
g!(zeros(3,1),Guess)
ForwardDiff.gradient(f,Guess)

h!(zeros(3,3),Guess)
ForwardDiff.hessian(f,Guess)

# ---------------------------#
# Exploring Optim #

# (1) Optimization using ForwardDiff: (working well)
Optim.minimizer(optimize(f, Guess)) # Same result as: (x'x)\(x'y)

# (2) Optimization using analytical gradient: (fail)
Optim.minimizer(optimize(f, g!, Guess)) # Approach 1

td = TwiceDifferentiable( f, g!, Guess; autodiff=:forward) # Approach 2
Optim.minimizer(optimize(td, Guess, LBFGS()))

# (2) Optimization using analytical hessian: (Issue: gives the same initial guess!)
Optim.minimizer(optimize(f, g!, h!, Guess, Newton()))

```

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [February 21, 2022, 7:51pm UTC](https://discourse.julialang.org/t/optim-user-supplied-gradient-and-hessian-mle/76848/2 "2022-02-21T19:51:00Z")

</div>

If I do

```julia
# prelim check that grad and hess are ok:
@show g!(zeros(3,1),Guess)
@show ForwardDiff.gradient(f,Guess)

@show h!(zeros(3,3),Guess)
@show ForwardDiff.hessian(f,Guess)

```

I see

```julia
g!(zeros(3, 1), Guess) = [1330.8310178125712; 1304.0590964682765; 4429.092663725574;;]
ForwardDiff.gradient(f, Guess) = [-1330.8310178125712, -1304.059096468278, -4429.092663725576]
h!(zeros(3, 3), Guess) = [-500.0 18.65451374962747 -1330.8310178125712; 18.65451374962747 -496.5116016197911 -1304.0590964682765; -1330.8310178125712 -1304.0590964682765 -9108.185327451149]
ForwardDiff.hessian(f, Guess) = [500.0 -18.654513749627476 1330.8310178125712; -18.654513749627476 496.5116016197912 1304.059096468278; 1330.8310178125712 1304.059096468278 9108.185327451152]

```

which looks slightly off? Furthermore `Optim.minimizer(optimize(f, g!, Guess))` fails with

```julia
ERROR: LoadError: DomainError with -4428.092663725574:

```

which can be directly generated with

```julia
g!(zeros(3), [-1329.8310178125712, -1303.0590964682765, -4428.092663725574])

```

---

<div class="post-metadata">

**Author:** ![luischanci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luischanci/32/29503_2.png) [@luischanci](https://discourse.julialang.org/u/luischanci)\
**Post date:** [February 22, 2022, 2:34am UTC](https://discourse.julialang.org/t/optim-user-supplied-gradient-and-hessian-mle/76848/3 "2022-02-22T02:34:51Z")

</div>

Thanks for pointing that out.  
The switch in the sign is related to the return of -LL, rather than LL, in the objective function. So, changing the sign in the returns of the functions for the analytical gradient and hessian may be enough to produce comparable results with those using ForwardDiff.  
On the other hand, your last point does suggest issues related to the specification of the function g!() - or, specifically, in ▽f. I will check that! (perhaps using exp(sigma2) rather than sigma2).  
In any case, I would like to explore the use of the hessian with Optim for the case of a normal distribution in MLE, so I would appreciate any advice.

---

<div class="post-metadata">

**Author:** ![luischanci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luischanci/32/29503_2.png) [@luischanci](https://discourse.julialang.org/u/luischanci)\
**Post date:** [February 22, 2022, 5:19pm UTC](https://discourse.julialang.org/t/optim-user-supplied-gradient-and-hessian-mle/76848/4 "2022-02-22T17:19:10Z")

</div>

To fix

```julia
DomainError with -4428.092663725574:
sqrt will only return a complex result if called with a complex argument.

```

I updated the functions to:

```julia
function Loglik(N,Y,X,Γ)
    σ² = exp(Γ[end])
    ep = (Y - X*Γ[1:(length(Γ)-1)])
    LL = -(N/2)*log(2π) - (N/2)*log(σ²) - (ep'*ep)/(2*σ²)
    return -LL
end
function ▽f(N,Y,X,▽,Γ)
    σ² = exp(Γ[end,1])
    k = (length(Γ)-1)
    ep = (Y - X*Γ[1:k])
    ▽[1:k] = (1/σ²)*(X'*ep)
    ▽[end] = - N/2 + (ep'*ep)/(2*σ²)
    return -▽
end
function ▽²f(N,Y,X,H,Γ)
    σ² = exp(Γ[end,1])
    k = (length(Γ)-1)
    ep = (Y - X*Γ[1:k])
    H[1:k,1:k] = -(1/σ²)*(X'X);
    H[1:k,end] = -(1/σ²)*(X'ep);
    H[end,1:k] = H[1:k,end]';
    H[end,end] = - (1/(2*σ²))*(ep'*ep)
    return -H
end

f(α) = Loglik(n,y,x,α)
g!(G,α) = ▽f(n,y,x,G,α)
h!(H,α) = ▽²f(n,y,x,H,α)
Guess = ones(size(x)[2]+1);

# prelim check that grad and hess are ok:
g!(zeros(3,1),Guess)
ForwardDiff.gradient(f,Guess)

h!(zeros(3,3),Guess)
ForwardDiff.hessian(f,Guess)

```

but if I do

```julia
Optim.minimizer(optimize(f, g!, Guess, LBFGS()))

```

I see the same initial guess

```julia
3-element Vector{Float64}:
 0.9995485513918926
 0.9995576333640301
 0.9986430247618775

```

likewise, if I do

```julia
Optim.minimizer(optimize(f, g!, h!, Guess, Newton()))

```

I see

```julia
3-element Vector{Float64}:
 0.9963506666680584
 0.996394580539055
 1.0007370702641363

```