The Hessian after optimization is not even close to symmetric; for example the 1, 2 element is 16.7 while the 2, 1 element is 1.4. The 12, 13 element is -3945 and 13, 12 is -16 (I’m rounding the numbers).
I’m using autodifferentiation, and so exactly where things are going wrong is unclear. These gross asymmetries appear only when
Newton() is the optimization method;
BFGS() does not (by eye) have them. Both methods produce essentially the same solution, though
BFGS reports convergence failure.
BFGS produces generally larger standard error estimates, and no negative variances.
num_hess in the code below is just a reference to
f, the latter gets garbage-collected, and a reference to random memory is left?
I’m using development code for StatsModels 0.6, though I haven’t updated it in a couple of weeks.
This may be related to my older, though still open, report of asymmetry (https://github.com/JuliaMath/Calculus.jl/issues/91), but that was in a much different context. And in my old report the asymmetry was small. See also https://github.com/JuliaDiff/ForwardDiff.jl/issues/253#issue-252103867.
This is not a nice, self-contained example, but the key code is
using DataFrames, StatsModels using NamedArrays using Optim, NLSolversBase using LinearAlgebra # for dot function mysolve(formula, df, weight=nothing)::AbstractResult nsamp = size(df, 1) if isnothing(weight) w = ones(nsamp) else w = nsamp/sum(weight)*weight end xp = prep(formula, df, w) myll(β) = -myllike(β, xp) f = TwiceDifferentiable(myll, zeros(size(xp.mm, 2)); autodiff=:forward) r = optimize(f, zeros(size(xp.mm, 2)), Newton(), Optim.Options(store_trace=true)) p = Optim.minimizer(r) num_hess = hessian!(f, p) # this is the Hessian I reported vcv = inv(num_hess) v = diag(vcv) sd = v sd[sd .< 0.0] .= 0.0 # hack for something that should not happen sd = sqrt.(sd) withname(x) = NamedArray(x, (xp.nms,)) with2name(x) = NamedArray(x, [xp.nms, xp.nms]) return SimpleResult(withname(p), withname(sd), with2name(vcv), with2name(num_hess), r, xp) end
More package versions:
[324d7699] CategoricalArrays v0.5.4 [a93c6f00] DataFrames v0.18.3 [f6369f11] ForwardDiff v0.10.3 [38e38edf] GLM v1.1.2 [`dev\GLM`] [54eb57ff] InteractiveCodeSearch v0.3.1 [d41bc354] NLSolversBase v7.3.1 [86f7a689] NamedArrays v0.9.3 [429524aa] Optim v0.18.1 [df47a6cb] RData v0.6.0 [3eaba693] StatsModels v0.6.0 [`dev\StatsModels`] [bd369af6] Tables v0.2.6
P.S. When I copy and paste out of the Junu REPL many line breaks vanish. Is there a way to fix that? I’m on Windows.