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.:

```
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.

```
## 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) )
```