Connection between BFGS and ROOT.Minuit. Stopping criteria

Specific Domains Optimization (Mathematical)
,
4

The stopping criteria can be implemented as follows.

"""
See also MIGrad in Chapter 4: Minuit Commands, https://root.cern.ch/download/minuit.pdf
"""
function miniutestop(state)
    mt = state.metadata
    edm = mt["g(x)"]' * mt["~inv(H)"] * mt["g(x)"] / 2
    edm < 1e-3 * 0.1 * 1.0
end

Here is a minimal working example:

using Distributions, Optim

"""
See also MIGrad in Chapter 4: Minuit Commands, https://root.cern.ch/download/minuit.pdf
""" 
function miunitstop(state)
    mt = state.metadata
    edm = mt["g(x)"]' * mt["~inv(H)"] * mt["g(x)"] / 2
    edm < 1e-3 * 0.1 * 1.0
end

data = randn(1000)
res = optimize([-Inf, 0.0], [Inf, Inf], [0.0, 1.0], Fminbox(BFGS()), Optim.Options(extended_trace=true, callback=miunitstop)) do pars
    -loglikelihood(Normal(pars...), data)
end

Please note that the callback has to be added in Optim.Options along with extended_trace=true for those metadata to be accessible.

Note also that the inverted hessian used in the calculation of EDM is just an approximate one that happens to be used in BFGS. So, after the minimum is found, you may want to calculate an exact one with Zygote.hessian(pars -> -loglikelihood(Normal(pars...), data), res.minimizer) (especially when you need to invert the hessian matrix to get the covariance matrix).

Demo

The miunitstop can reduce the number of calls by ~100x, with effectively the same results:

julia> res = optimize([-Inf, 0.0], [Inf, Inf], [0.0, 1.0], Fminbox(BFGS()), Optim.Options(extended_trace=true, callback=miunitstop)) do pars
           -loglikelihood(Normal(pars...), data)
       end
 * Status: failure

 * Candidate solution
    Final objective value:     1.421423e+03

 * Found with
    Algorithm:     Fminbox with BFGS

 * Convergence measures
    |x - x'|               = 1.34e-02 ≰ 0.0e+00
    |x - x'|/|x'|          = 1.34e-02 ≰ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 2.06e-02 ≰ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    1
    f(x) calls:    9
    ∇f(x) calls:   9


julia> res2 = optimize([-Inf, 0.0], [Inf, Inf], [0.0, 1.0], Fminbox(BFGS())) do pars
           -loglikelihood(Normal(pars...), data)
       end
 * Status: success (objective increased between iterations)

 * Candidate solution
    Final objective value:     1.421423e+03

 * Found with
    Algorithm:     Fminbox with BFGS

 * Convergence measures
    |x - x'|               = 9.17e-11 ≰ 0.0e+00
    |x - x'|/|x'|          = 9.15e-11 ≰ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 3.75e-08 ≰ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    4
    f(x) calls:    713
    ∇f(x) calls:   713

julia> res.minimizer
2-element Vector{Float64}:
 0.013204458490624406
 1.002498042266284

julia> res2.minimizer
2-element Vector{Float64}:
 0.013199480602918739
 1.002487701716321
5 Likes