Wrongly reporting success

[azev77]

I'm solving for the MLE.

using Distributions, Random, Optim;
DGP_True = LogNormal(17,7);
Random.seed!(123);
const d_train = rand(DGP_True, 1_000)
const d_test  = rand(DGP_True, 1_000)

# LogNormal
function loglike(θ)
    dist          = LogNormal(θ[1],θ[2])
    terms         = logpdf.(dist,d_train)
    loglikelihood = sum(terms)
    return -loglikelihood #Min(-LogLik) =Max(LogLik)
end
#
optimum = optimize(loglike,ones(2))
Optim.minimizer(optimum)
Distributions.fit_mle(LogNormal,d_train)

# Optim.jl &  Distributions.jl solutions are very close.


# Normal 
function loglike(θ)
    dist          = Normal(θ[1],θ[2])
    terms         = logpdf.(dist,d_train)
    loglikelihood = sum(terms)
    return -loglikelihood #Min(-LogLik) =Max(LogLik)
end
#
optimum = optimize(loglike,ones(2))
Optim.minimizer(optimum)
Distributions.fit_mle(Normal,d_train)

# Optim.jl reports "sucess" BUT
# Optim.jl &  Distributions.jl solutions are very different!

[azev77]

Here is a 2nd bug:

using Distributions, Random, Optim;
DGP_True = LogNormal(0,1);
Random.seed!(123);
const d_train = rand(DGP_True, 1_000);
const d_test  = rand(DGP_True, 1_000);

D1=Distributions.fit_mle(Laplace, d_train)

function loglike(θ)
    dist          = Laplace(θ[1],θ[2])
    terms         = logpdf.(dist,d_train)
    loglikelihood = sum(terms)
    return -loglikelihood #Min(-LogLik) =Max(LogLik)
end
sol=optimize(loglike,ones(2),Optim.Options(iterations=10^4))
D2=Optim.minimizer(sol)

The estimate for the first parameter is very similar:
μ=1.0112216426398657 vs Optim's 1.0116066436206388
The estimate for the 2nd parameter is different:
θ=0.9005083163386765 vs Optim's 1.1872092883101257

This is concerning bc Optim is reporting a "success".

[johnmyleswhite]

This is an interesting example and hopefully Optim can do more here. The problem seems to be that your log likelihood function essentially behaves like a constant over most reasonable regions of parameter space and Optim doesn't learn quickly enough that your problem is so poorly scaled and doesn't know how to compensate for the scaling problems.

Since your problem is differentiable, it's usually helpful to use L-BFGS; if you do that, you quickly run into an issue where the gradient after a single step of minimization is essentially zero (and it's even smaller than the gradient at the analytic MLE). If you tell Optim to ignore that gradient information, you then immediately hit a line search failure:

julia> using Distributions, Random, Optim;

julia> DGP_True = LogNormal(17,7);

julia> Random.seed!(123);

julia> const d_train = rand(DGP_True, 1_000);

julia> const d_test  = rand(DGP_True, 1_000);

julia> import ForwardDiff: gradient

julia> function loglike(θ)
           dist          = Normal(θ[1],θ[2])
           terms         = logpdf.(dist,d_train)
           loglikelihood = sum(terms)
           return -loglikelihood #Min(-LogLik) =Max(LogLik)
       end
loglike (generic function with 1 method)

julia> optimum = optimize(
           loglike,
           ones(2),
           Optim.LBFGS(),
           autodiff = :forward,
           Optim.Options(g_tol = 0.0, show_trace = true)
       )
Iter     Function value   Gradient norm
     0     1.061972e+32     2.123944e+32
 * time: 9.107589721679688e-5
     1     7.535494e+04     4.708222e-30
 * time: 0.7484478950500488
 * Status: failure (line search failed)

 * Candidate solution
    Final objective value:     7.535494e+04

 * Found with
    Algorithm:     L-BFGS

 * Convergence measures
    |x - x'|               = 0.00e+00 ≤ 0.0e+00
    |x - x'|/|x'|          = 0.00e+00 ≤ 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)|                 = 4.71e-30 ≰ 0.0e+00

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

julia> (
           loglike(Optim.minimizer(optimum))
           ===
           loglike(Optim.minimizer(optimum) .+ eps(Optim.minimizer(optimum)[1]))
       )
true

julia> (
           loglike(Optim.minimizer(optimum))
           ===
           loglike(Optim.minimizer(optimum) .+ 1000000 * eps(Optim.minimizer(optimum)[1]))
       )
true

julia> import ForwardDiff: gradient

julia> gradient(loglike, Optim.minimizer(optimum))
2-element Array{Float64,1}:
 8.127227357617317e-46
 4.708222339345497e-30

julia> gradient(loglike, [mean(d_train), std(d_train)])
2-element Array{Float64,1}:
 -8.578862344278503e-30
  2.1756679847282367e-15

It's not totally obvious to me how Optim can work around the log likelihood function being problematic numerically.

[timholy]

There are very few convergence criteria that pass the sniff test of being mathematically reasonable. So many of them can be defeated simply by defining

fsc(args...) = 1e-30*f(args...)

and then running the optimization on fsc. But obviously multiplying by a constant should have no impact on the definition of the location of the minimum.

If we didn't have to worry about objectives having a minimum value of exactly 0.0, this would not be a problem; we could use all relative criteria, i.e., abs(fold-fnew) <= tol*(abs(fold) + abs(fnew)). But given that some functions do have a minimum value of 0, there doesn't seem to be an easy solution.

[antoine-levitt]

The solution would be to pass a reasonable scale of your problem to Optim (both for f and for x). But then if you know a reasonable scale you probably want to scale f and x yourself, so I don't think it would make much sense for Optim to do this. Is the problem in the OP simply a question of scaling, or is it something like x -> exp(-(x-x*)^2 / 100000)? If the latter, this is a very challenging optimization problem that can't really be solved by a uniform scaling...

[pkofod]

This is concerning bc Optim is reporting a "success".

Yes, I agree, I don't like that it ever prints success. It should just print the different convergence information.

I'm wondering what the difference between the trace might be if you switch to automatic differentiation like John did. But in the end there's only so much you can do with poor scaling. If it's not completely constant just poorly scaled you can try to use Newton's method.

[johnmyleswhite]

AD just exposes the issue that you can't have "nearly" zero gradient check as a convergence criteria if the function is constant in a small region that it might reach. You can just remove the explicit g_tol = 0.0 and see this play out:

julia> optimum = optimize(
           loglike,
           ones(2),
           Optim.LBFGS(),
           autodiff = :forward,
           Optim.Options(show_trace = true)
       )
Iter     Function value   Gradient norm
     0     1.061972e+32     2.123944e+32
 * time: 7.295608520507812e-5
     1     7.535494e+04     4.708222e-30
 * time: 0.5777468681335449
 * Status: success

 * Candidate solution
    Final objective value:     7.535494e+04

 * Found with
    Algorithm:     L-BFGS

 * Convergence measures
    |x - x'|               = 2.12e+32 ≰ 0.0e+00
    |x - x'|/|x'|          = 1.00e+00 ≰ 0.0e+00
    |f(x) - f(x')|         = 1.06e+32 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 1.41e+27 ≰ 0.0e+00
    |g(x)|                 = 4.71e-30 ≤ 1.0e-08

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

I'm really not clear what Optim can do here -- the function is just really badly behaved.

Maybe it's useful to ask what Optim would do for the function defined:

• f(x) = 1 if x <= -1
• f(x) = x^2 if -1 <= x <= 1
• f(x) = 1 if 1 <= x

If you start that at, say, x = 10, you're going to get success because you are at a local minimizer, no?

julia> optimize(f, [-10.0], LBFGS())
 * Status: success

 * Candidate solution
    Final objective value:     1.000000e+00

 * Found with
    Algorithm:     L-BFGS

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

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

One could replace "success" with "converged = true", but I'm not sure that's enough to fix the issue that some users will send Optim a function that isn't really tractable to the convergence criteria in use.