Box-constrained optimization calls values outside of the box

I’m trying to use Optim.jl to minimize a function that has hard boundaries on the input using the Fminbox box-constrained optimization. When the minimum value is near or on the boundary, however, the optimizer can check values outside the defined box. I’m assuming it’s because the boundaries aren’t hard constraints, but rather there’s a barrier function that’s insufficient for keeping the optimizer inside the box. I don’t know how to overcome that though.

Here’s a MWE in 1D. This problem could be solved with a univariate optimizer, of course, but my actual application is multivariate:

using Optim
function boxmin(xmin)
    x0 = [0.1]
    lower = [0.0]
    upper = [1.0]
    function f(x)
        @assert x[1] >= lower[1]
        @assert x[1] <= upper[1]
        return (x[1]-xmin)^2
    end
    res = Optim.optimize(f, lower, upper, x0, Optim.Fminbox(Optim.GradientDescent()),Optim.Options(store_trace = true, show_trace = true))
end

where xmin is the location of the mimimum.

For xmin=1e-4, this will converge:

Fminbox
-------
Initial mu = 1.7982e-5

Fminbox iteration 1
-------------------
Calling inner optimizer with mu = 1.7982e-5

(numbers below include barrier contribution)
Iter     Function value   Gradient norm 
     0     1.002331e-02     1.998433e-01
 * time: 0.00011491775512695312
     1     4.339978e-04     3.715841e-02
 * time: 0.0006890296936035156
     2     1.137071e-04     2.346214e-03
 * time: 0.0011658668518066406
     3     1.129297e-04     2.724118e-04
 * time: 0.0015418529510498047
     4     1.129201e-04     2.024939e-06
 * time: 0.0019528865814208984
     5     1.129201e-04     1.664494e-11
 * time: 0.002346038818359375

Exiting inner optimizer with x = [0.0030443354558911933]
Current distance to box: 0.00304434
Decreasing barrier term μ.

Fminbox iteration 2
-------------------
Calling inner optimizer with mu = 1.7982e-8

(numbers below include barrier contribution)
Iter     Function value   Gradient norm 
     0     8.773362e-06     5.888775e-03
 * time: 3.2901763916015625e-5
     1     1.648191e-07     1.419758e-04
 * time: 0.00042700767517089844
     2     1.607654e-07     4.950859e-06
 * time: 0.0007879734039306641
     3     1.607609e-07     7.838630e-09
 * time: 0.0011649131774902344

Exiting inner optimizer with x = [0.00015719236235675393]
Current distance to box: 0.000157192
Decreasing barrier term μ.

Fminbox iteration 3
-------------------
Calling inner optimizer with mu = 1.7982e-11

(numbers below include barrier contribution)
Iter     Function value   Gradient norm 
     0     3.428456e-09     1.143849e-04
 * time: 2.8848648071289062e-5
     1     1.656141e-10     1.118057e-08
 * time: 0.0004899501800537109
     2     1.656141e-10     1.877634e-18
 * time: 0.00096893310546875

Exiting inner optimizer with x = [0.00010008982032311996]
Current distance to box: 0.00010009
Decreasing barrier term μ.

Fminbox iteration 4
-------------------
Calling inner optimizer with mu = 1.7982e-14

(numbers below include barrier contribution)
Iter     Function value   Gradient norm 
     0     1.736737e-13     1.796408e-07
 * time: 4.696846008300781e-5
     1     1.656221e-13     1.239897e-19
 * time: 0.0005347728729248047

Exiting inner optimizer with x = [0.00010000008990092734]
Current distance to box: 0.0001
Decreasing barrier term μ.

 * Status: success

 * Candidate solution
    Final objective value:     8.082177e-21

 * Found with
    Algorithm:     Fminbox with Gradient Descent

 * Convergence measures
    |x - x'|               = 8.97e-08 ≰ 0.0e+00
    |x - x'|/|x'|          = 8.97e-04 ≰ 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)|                 = 1.80e-10 ≤ 1.0e-08

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

For xmin=1e-5, this fails:

Fminbox
-------
Initial mu = 1.79982e-5

Fminbox iteration 1
-------------------
Calling inner optimizer with mu = 1.79982e-5

(numbers below include barrier contribution)
Iter     Function value   Gradient norm 
     0     1.004134e-02     2.000233e-01
 * time: 8.797645568847656e-5
     1     6.200988e-04     4.624942e-02
 * time: 0.0005919933319091797
     2     1.139175e-04     1.645394e-03
 * time: 0.0011169910430908203
     3     1.135500e-04     1.179194e-04
 * time: 0.0015099048614501953
     4     1.135483e-04     2.980599e-07
 * time: 0.0018749237060546875
     5     1.135483e-04     1.548770e-13
 * time: 0.0021889209747314453

Exiting inner optimizer with x = [0.00300033694501183]
Current distance to box: 0.00300034
Decreasing barrier term μ.

Fminbox iteration 2
-------------------
Calling inner optimizer with mu = 1.79982e-8

(numbers below include barrier contribution)
Iter     Function value   Gradient norm 
     0     9.046721e-06     5.980778e-03
 * time: 3.1948089599609375e-5
     1     1.742513e-07     5.160145e-05
 * time: 0.00037384033203125
     2     1.738728e-07     3.305504e-06
 * time: 0.0007488727569580078
     3     1.738713e-07     6.676051e-09
 * time: 0.001071929931640625

Exiting inner optimizer with x = [9.999228308793515e-5]
Current distance to box: 9.99923e-5
Decreasing barrier term μ.

Fminbox iteration 3
-------------------
Calling inner optimizer with mu = 1.79982e-11

(numbers below include barrier contribution)
Iter     Function value   Gradient norm 
     0     8.264384e-09     1.799847e-04
 * time: 6.29425048828125e-5
AssertionError: x[1] >= lower[1]

Stacktrace:
  [1] f
    @ ./In[18]:9 [inlined]
  [2] finite_difference_gradient!(df::Vector{Float64}, f::var"#f#4"{Float64, Vector{Float64}, Vector{Float64}}, x::Vector{Float64}, cache::FiniteDiff.GradientCache{Nothing, Nothing, Nothing, Vector{Float64}, Val{:central}(), Float64, Val{true}()}; relstep::Float64, absstep::Float64, dir::Bool)
    @ FiniteDiff ~/.julia/packages/FiniteDiff/40JnL/src/gradients.jl:320
  [3] finite_difference_gradient!
    @ ~/.julia/packages/FiniteDiff/40JnL/src/gradients.jl:258 [inlined]
  [4] (::NLSolversBase.var"#g!#15"{var"#f#4"{Float64, Vector{Float64}, Vector{Float64}}, FiniteDiff.GradientCache{Nothing, Nothing, Nothing, Vector{Float64}, Val{:central}(), Float64, Val{true}()}})(storage::Vector{Float64}, x::Vector{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/objective_types/oncedifferentiable.jl:57
  [5] (::NLSolversBase.var"#fg!#16"{var"#f#4"{Float64, Vector{Float64}, Vector{Float64}}})(storage::Vector{Float64}, x::Vector{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/objective_types/oncedifferentiable.jl:61
  [6] value_gradient!!(obj::OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, x::Vector{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/interface.jl:82
  [7] value_gradient!(obj::OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, x::Vector{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/interface.jl:69
  [8] value_gradient!(bb::Optim.BarrierWrapper{OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BoxBarrier{Vector{Float64}, Vector{Float64}}, Float64, Float64, Vector{Float64}}, x::Vector{Float64})
    @ Optim ~/.julia/packages/Optim/tP8PJ/src/multivariate/solvers/constrained/fminbox.jl:94
  [9] value_gradient!(obj::Optim.ManifoldObjective{Optim.BarrierWrapper{OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BoxBarrier{Vector{Float64}, Vector{Float64}}, Float64, Float64, Vector{Float64}}}, x::Vector{Float64})
    @ Optim ~/.julia/packages/Optim/tP8PJ/src/Manifolds.jl:50
 [10] (::LineSearches.var"#ϕdϕ#6"{Optim.ManifoldObjective{Optim.BarrierWrapper{OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BoxBarrier{Vector{Float64}, Vector{Float64}}, Float64, Float64, Vector{Float64}}}, Vector{Float64}, Vector{Float64}, Vector{Float64}})(α::Float64)
    @ LineSearches ~/.julia/packages/LineSearches/G1LRk/src/LineSearches.jl:84
 [11] (::LineSearches.HagerZhang{Float64, Base.RefValue{Bool}})(ϕ::Function, ϕdϕ::LineSearches.var"#ϕdϕ#6"{Optim.ManifoldObjective{Optim.BarrierWrapper{OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BoxBarrier{Vector{Float64}, Vector{Float64}}, Float64, Float64, Vector{Float64}}}, Vector{Float64}, Vector{Float64}, Vector{Float64}}, c::Float64, phi_0::Float64, dphi_0::Float64)
    @ LineSearches ~/.julia/packages/LineSearches/G1LRk/src/hagerzhang.jl:235
 [12] HagerZhang
    @ ~/.julia/packages/LineSearches/G1LRk/src/hagerzhang.jl:101 [inlined]
 [13] perform_linesearch!(state::Optim.GradientDescentState{Vector{Float64}, Float64}, method::GradientDescent{LineSearches.InitialPrevious{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Optim.InverseDiagonal, Optim.var"#66#67"{Vector{Float64}, Vector{Float64}, Fminbox{GradientDescent{LineSearches.InitialPrevious{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Optim.var"#13#15"}, Float64, Optim.var"#49#51"}, Optim.BarrierWrapper{OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BoxBarrier{Vector{Float64}, Vector{Float64}}, Float64, Float64, Vector{Float64}}}}, d::Optim.ManifoldObjective{Optim.BarrierWrapper{OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BoxBarrier{Vector{Float64}, Vector{Float64}}, Float64, Float64, Vector{Float64}}})
    @ Optim ~/.julia/packages/Optim/tP8PJ/src/utilities/perform_linesearch.jl:59
 [14] update_state!(d::Optim.BarrierWrapper{OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BoxBarrier{Vector{Float64}, Vector{Float64}}, Float64, Float64, Vector{Float64}}, state::Optim.GradientDescentState{Vector{Float64}, Float64}, method::GradientDescent{LineSearches.InitialPrevious{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Optim.InverseDiagonal, Optim.var"#66#67"{Vector{Float64}, Vector{Float64}, Fminbox{GradientDescent{LineSearches.InitialPrevious{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Optim.var"#13#15"}, Float64, Optim.var"#49#51"}, Optim.BarrierWrapper{OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BoxBarrier{Vector{Float64}, Vector{Float64}}, Float64, Float64, Vector{Float64}}}})
    @ Optim ~/.julia/packages/Optim/tP8PJ/src/multivariate/solvers/first_order/gradient_descent.jl:80
 [15] optimize(d::Optim.BarrierWrapper{OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BoxBarrier{Vector{Float64}, Vector{Float64}}, Float64, Float64, Vector{Float64}}, initial_x::Vector{Float64}, method::GradientDescent{LineSearches.InitialPrevious{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Optim.InverseDiagonal, Optim.var"#66#67"{Vector{Float64}, Vector{Float64}, Fminbox{GradientDescent{LineSearches.InitialPrevious{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Optim.var"#13#15"}, Float64, Optim.var"#49#51"}, Optim.BarrierWrapper{OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BoxBarrier{Vector{Float64}, Vector{Float64}}, Float64, Float64, Vector{Float64}}}}, options::Optim.Options{Float64, Nothing}, state::Optim.GradientDescentState{Vector{Float64}, Float64})
    @ Optim ~/.julia/packages/Optim/tP8PJ/src/multivariate/optimize/optimize.jl:54
 [16] optimize(df::OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, l::Vector{Float64}, u::Vector{Float64}, initial_x::Vector{Float64}, F::Fminbox{GradientDescent{LineSearches.InitialPrevious{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Optim.var"#13#15"}, Float64, Optim.var"#49#51"}, options::Optim.Options{Float64, Nothing})
    @ Optim ~/.julia/packages/Optim/tP8PJ/src/multivariate/solvers/constrained/fminbox.jl:390
 [17] #optimize#65
    @ ~/.julia/packages/Optim/tP8PJ/src/multivariate/solvers/constrained/fminbox.jl:269 [inlined]
 [18] optimize
    @ ~/.julia/packages/Optim/tP8PJ/src/multivariate/solvers/constrained/fminbox.jl:259 [inlined]
 [19] boxmin(xmin::Float64)
    @ Main ./In[18]:14
 [20] top-level scope
    @ In[23]:1
2 Likes

I’m not familiar with the Optim.jl code, but this suggests that the functions are evaluated before the fraction-to-boundary rule (responsible for enforcing strict bound constraints) is enforced. That’s weird… I’ll check the code.