# Why is Optim.optimize() throwing an exception if I provide initial\_invH?

**URL:** <https://discourse.julialang.org/t/why-is-optim-optimize-throwing-an-exception-if-i-provide-initial-invh/59872>\
**Category:** Optimization (Mathematical)\
**Tags:** question\
**Created:** [April 23, 2021, 12:41pm UTC](https://discourse.julialang.org/t/why-is-optim-optimize-throwing-an-exception-if-i-provide-initial-invh/59872 "2021-04-23T12:41:37Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![marcus-waldman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marcus-waldman/32/24196_2.png) [@marcus-waldman](https://discourse.julialang.org/u/marcus-waldman)\
**Post date:** [April 23, 2021, 12:41pm UTC](https://discourse.julialang.org/t/why-is-optim-optimize-throwing-an-exception-if-i-provide-initial-invh/59872/1 "2021-04-23T12:41:37Z")

</div>

I’m trying to implement the EM algorithm to estimate parameters in a large model The EM algorithm requires a lot of back and forth updating some parameters of interest and then optimizing an objective function based on the current state of the parameter values. The latter step I am completing using `Optim.optimize()`.

I am trying to provide a “current” best guess of the (inverse) Hessian matrix when calling `Optim.optimize()` and using `BFGS()`. The specific problem I have is that passing this best-guess initial Hessian results in a methods exception being thrown, i.e.:

```julia
opt2 = optimize(func, Optim.minimizer(opt1), 
                BFGS(initial_invH = opt1.trace[end].metadata["~inv(H)"]), 
                Optim.Options(show_trace = true) )

#note opt1 comes from the optimization at a previous iteration in the routine. 

```

results in

> MethodError: objects of type Matrix{Float64} are not callable  
> Use square brackets for indexing an Array.  
> initial\_state(method::BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Matrix{Float64}, Nothing, Flat}, options::Optim.Options{Float64, Nothing}, d::OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, initial\_x::Vector{Float64}) at bfgs.jl:97  
> optimize(d::OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, initial\_x::Vector{Float64}, method::BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Matrix{Float64}, Nothing, Flat}, options::Optim.Options{Float64, Nothing}) at optimize.jl:35  
> top-level scope at mwe-provide-invH.jl:28  
> eval at boot.jl:360 [inlined]

A minimum working example is below for reproducibility. It is not the EM algorithm I am implementing in my actual code, but the same exception pops up. Any assistance is, of course, greatly appreciated as being able to provide an initial Hessian would greatly speed up the execution of my code.

```julia
## Following minimum working example is a modified version of what appears in 
# https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/maxlikenlm/
using Optim, NLSolversBase, Random
using LinearAlgebra: diag
Random.seed!(0)

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

function Log_Likelihood(X, Y, β, log_σ)
    σ = exp(log_σ)
    llike = -n/2*log(2π) - n/2* log(σ^2) - (sum((Y - X * β).^2) / (2σ^2))
    llike = -llike
end

func = OnceDifferentiable(vars -> Log_Likelihood(x, y, vars[1:nvar], vars[nvar + 1]),
                           ones(nvar+1); autodiff=:forward);
                           
opt1 = optimize(func, ones(nvar+1), 
                BFGS(), 
                Optim.Options(store_trace = true, extended_trace=true, iterations = 3))

opt2 = optimize(func, Optim.minimizer(opt1), 
                BFGS(initial_invH = opt1.trace[end].metadata["~inv(H)"]), 
                Optim.Options(show_trace = true) )

```

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 23, 2021, 1:40pm UTC](https://discourse.julialang.org/t/why-is-optim-optimize-throwing-an-exception-if-i-provide-initial-invh/59872/2 "2021-04-23T13:40:31Z")

</div>

Thanks for taking the time to construct an MWE. I think that constructor expects a closure (you can provide an initial inverse Hessian that depends on x), so something like

```julia
opt2 = optimize(func, Optim.minimizer(opt1),
                BFGS(initial_invH = _ -> opt1.trace[end].metadata["~inv(H)"]),
                Optim.Options(show_trace = true) )

```

would work (note the `_ ->`).

---

<div class="post-metadata">

**Author:** ![marcus-waldman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marcus-waldman/32/24196_2.png) [@marcus-waldman](https://discourse.julialang.org/u/marcus-waldman)\
**Post date:** [April 23, 2021, 2:30pm UTC](https://discourse.julialang.org/t/why-is-optim-optimize-throwing-an-exception-if-i-provide-initial-invh/59872/3 "2021-04-23T14:30:23Z")

</div>

Thanks @Tamas_Papp. This worked. Now that I know what I am looking for, I can see how you arrived at that from the documentation `?BFGS()`.

> # Constructor
> 
> BFGS(; alphaguess = LineSearches.InitialStatic(),  
> linesearch = LineSearches.HagerZhang(),  
> initial\_invH = x → Matrix{eltype(x)}(I, length(x), length(x)),  
> manifold = Flat())

I wonder if there is a way for applied users to contribute to the documentation to provide concrete examples within the documentation, so as to elaborate on arguments. I really like Julia and would love to promote it to applied researchers in my field as an alternative to R, but I think they would have a difficult time without more explicit handholding in the documentation and help files. I know developers are tied up on other important things, however.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 25, 2021, 8:08am UTC](https://discourse.julialang.org/t/why-is-optim-optimize-throwing-an-exception-if-i-provide-initial-invh/59872/4 "2021-04-25T08:08:27Z")

</div>

> [@marcus-waldman](#):
>
> I wonder if there is a way for applied users to contribute to the documentation to provide concrete examples within the documentation

Sure, just make a pull request on Github. If you only want to edit docstrings, maybe you can do that using the web interface.

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [April 25, 2021, 10:06am UTC](https://discourse.julialang.org/t/why-is-optim-optimize-throwing-an-exception-if-i-provide-initial-invh/59872/5 "2021-04-25T10:06:59Z")

</div>

The contribution would be very welcome 🙂