[WIP] Do not (mis)use objective as state

Author: devmotion

Project.toml

@@ -30,7 +30,7 @@ ExplicitImports = "1.13.2"
 FillArrays = "0.6.2, 0.7, 0.8, 0.9, 0.10, 0.11, 0.12, 0.13, 1"
 ForwardDiff = "0.10, 1"
 JET = "0.9, 0.10"
-LineSearches = "7.4.0"
+LineSearches = "7.5.1"
 LinearAlgebra = "<0.0.1, 1.6"
 MathOptInterface = "1.17"
 Measurements = "2.14.1"

src/Manifolds.jl

@@ -28,16 +28,7 @@ end
 # TODO: is it safe here to call retract! and change x?
 function NLSolversBase.value!(obj::ManifoldObjective, x)
     xin = retract(obj.manifold, x)
-    value!(obj.inner_obj, xin)
-end
-function NLSolversBase.value(obj::ManifoldObjective)
-    value(obj.inner_obj)
-end
-function NLSolversBase.gradient(obj::ManifoldObjective)
-    gradient(obj.inner_obj)
-end
-function NLSolversBase.gradient(obj::ManifoldObjective, i::Int)
-    gradient(obj.inner_obj, i)
+    return value!(obj.inner_obj, xin)
 end
 function NLSolversBase.gradient!(obj::ManifoldObjective, x)
     xin = retract(obj.manifold, x)

src/Optim.jl

@@ -41,10 +41,7 @@ using NLSolversBase:
     TwiceDifferentiableConstraints,
     nconstraints,
     nconstraints_x,
-    hessian,
     hessian!,
-    hessian!!,
-    hv_product,
     hv_product!
 
 # var for NelderMead

src/multivariate/optimize/optimize.jl

@@ -1,37 +1,57 @@
-update_g!(d, state, method) = nothing
+update_g!(d, state, ::ZerothOrderOptimizer) = nothing
 function update_g!(d, state, method::FirstOrderOptimizer)
     # Update the function value and gradient
-    value_gradient!(d, state.x)
-    project_tangent!(method.manifold, gradient(d), state.x)
+    f_x, g_x = value_gradient!(d, state.x)
+    project_tangent!(method.manifold, g_x, state.x)
+    state.f_x = f_x
+    copyto!(state.g_x, g_x)
+    return nothing
 end
-function update_g!(d, state, method::Newton)
+function update_g!(d, state, ::Newton)
     # Update the function value and gradient
-    value_gradient!(d, state.x)
+    f_x, g_x = value_gradient!(d, state.x)
+    state.f_x = f_x
+    copyto!(state.g_x, g_x)
+    return nothing
 end
-update_fg!(d, state, method) = nothing
-update_fg!(d, state, method::ZerothOrderOptimizer) = value!(d, state.x)
-function update_fg!(d, state, method::FirstOrderOptimizer)
-    value_gradient!(d, state.x)
-    project_tangent!(method.manifold, gradient(d), state.x)
+
+function update_fg!(d, state, method::ZerothOrderOptimizer)
+    f_x = value!(d, state.x)
+    state.f_x = f_x
+    return nothing
+end
+function update_fg!(d, state, method)
+    f_x, g_x = value_gradient!(d, state.x)
+    project_tangent!(method.manifold, g_x, state.x)
+    state.f_x = f_x
+    copyto!(state.g_x, g_x)
+    return nothing
 end
 function update_fg!(d, state, method::Newton)
-    value_gradient!(d, state.x)
+    f_x, g_x = value_gradient!(d, state.x)
+    state.f_x = f_x
+    copyto!(state.g_x, g_x)
+    return nothing
 end
 
 # Update the Hessian
-update_h!(d, state, method) = nothing
-update_h!(d, state, method::SecondOrderOptimizer) = hessian!(d, state.x)
+update_h!(d, state, method::Union{ZerothOrderOptimizer,FirstOrderOptimizer}) = nothing
+function update_h!(d, state, method::SecondOrderOptimizer)
+    H_x = hessian!(d, state.x)
+    copyto!(state.H_x, H_x)
+    return nothing
+end
 
 after_while!(d, state, method, options) = nothing
 
-function initial_convergence(d, state, method::AbstractOptimizer, initial_x, options)
-    gradient!(d, initial_x)
-    stopped = !isfinite(value(d)) || any(!isfinite, gradient(d))
-    g_residual(d, state) <= options.g_abstol, stopped
+function initial_convergence(state::AbstractOptimizerState, options::Options)
+    stopped = !isfinite(state.f_x) || any(!isfinite, state.g_x)
+    return g_residual(state) <= options.g_abstol, stopped
 end
-function initial_convergence(d, state, method::ZerothOrderOptimizer, initial_x, options)
+function initial_convergence(::ZerothOrderState, ::Options)
     false, false
 end
+
 function optimize(
     d::D,
     initial_x::Tx,
@@ -41,7 +61,7 @@ function optimize(
 ) where {D<:AbstractObjective,M<:AbstractOptimizer,Tx<:AbstractArray,T,TCallback}
 
     t0 = time() # Initial time stamp used to control early stopping by options.time_limit
-    tr = OptimizationTrace{typeof(value(d)),typeof(method)}()
+    tr = OptimizationTrace{typeof(state.f_x),typeof(method)}()
     tracing =
         options.store_trace ||
         options.show_trace ||
@@ -51,7 +71,7 @@ function optimize(
     f_limit_reached, g_limit_reached, h_limit_reached = false, false, false
     x_converged, f_converged, f_increased, counter_f_tol = false, false, false, 0
 
-    g_converged, stopped = initial_convergence(d, state, method, initial_x, options)
+    g_converged, stopped = initial_convergence(state, options)
     converged = g_converged || stopped
     # prepare iteration counter (used to make "initial state" trace entry)
     iteration = 0
@@ -113,11 +133,11 @@ function optimize(
             end
         end
 
-        if g_calls(d) > 0 && !all(isfinite, gradient(d))
+        if hasproperty(state, :g_x) && !all(isfinite, state.g_x)
             options.show_warnings && @warn "Terminated early due to NaN in gradient."
             break
         end
-        if h_calls(d) > 0 && !(d isa TwiceDifferentiableHV) && !all(isfinite, hessian(d))
+        if hasproperty(state, :H_x) && !all(isfinite, state.H_x)
             options.show_warnings && @warn "Terminated early due to NaN in Hessian."
             break
         end
@@ -127,7 +147,7 @@ function optimize(
 
     # we can just check minimum, as we've earlier enforced same types/eltypes
     # in variables besides the option settings
-    Tf = typeof(value(d))
+    Tf = typeof(state.f_x)
     f_incr_pick = f_increased && !options.allow_f_increases
     stopped_by = (x_converged, f_converged, g_converged,
         f_limit_reached = f_limit_reached,
@@ -141,7 +161,7 @@ function optimize(
     )
 
     termination_code =
-        _termination_code(d, g_residual(d, state), state, stopped_by, options)
+        _termination_code(d, g_residual(state), state, stopped_by, options)
 
     return MultivariateOptimizationResults{
         typeof(method),
@@ -154,18 +174,18 @@ function optimize(
         method,
         initial_x,
         pick_best_x(f_incr_pick, state),
-        pick_best_f(f_incr_pick, state, d),
+        pick_best_f(f_incr_pick, state),
         iteration,
         Tf(options.x_abstol),
         Tf(options.x_reltol),
         x_abschange(state),
         x_relchange(state),
         Tf(options.f_abstol),
         Tf(options.f_reltol),
-        f_abschange(d, state),
-        f_relchange(d, state),
+        f_abschange(state),
+        f_relchange(state),
         Tf(options.g_abstol),
-        g_residual(d, state),
+        g_residual(state),
         tr,
         f_calls(d),
         g_calls(d),
@@ -186,13 +206,13 @@ function _termination_code(d, gres, state, stopped_by, options)
     elseif (iszero(options.x_abstol) && x_abschange(state) <= options.x_abstol) ||
            (iszero(options.x_reltol) && x_relchange(state) <= options.x_reltol)
         TerminationCode.NoXChange
-    elseif (iszero(options.f_abstol) && f_abschange(d, state) <= options.f_abstol) ||
-           (iszero(options.f_reltol) && f_relchange(d, state) <= options.f_reltol)
+    elseif (iszero(options.f_abstol) && f_abschange(state) <= options.f_abstol) ||
+           (iszero(options.f_reltol) && f_relchange(state) <= options.f_reltol)
         TerminationCode.NoObjectiveChange
     elseif x_abschange(state) <= options.x_abstol || x_relchange(state) <= options.x_reltol
         TerminationCode.SmallXChange
-    elseif f_abschange(d, state) <= options.f_abstol ||
-           f_relchange(d, state) <= options.f_reltol
+    elseif f_abschange(state) <= options.f_abstol ||
+           f_relchange(state) <= options.f_reltol
         TerminationCode.SmallObjectiveChange
     elseif stopped_by.ls_failed
         TerminationCode.FailedLinesearch
@@ -210,11 +230,11 @@ function _termination_code(d, gres, state, stopped_by, options)
         TerminationCode.HessianCalls
     elseif stopped_by.f_increased
         TerminationCode.ObjectiveIncreased
-    elseif f_calls(d) > 0 && !isfinite(value(d))
-        TerminationCode.GradientNotFinite
-    elseif g_calls(d) > 0 && !all(isfinite, gradient(d))
+    elseif !isfinite(state.f_x)
+        TerminationCode.ObjectiveNotFinite
+    elseif hasproperty(state, :g_x) && !all(isfinite, state.g_x)
         TerminationCode.GradientNotFinite
-    elseif h_calls(d) > 0 && !(d isa TwiceDifferentiableHV) && !all(isfinite, hessian(d))
+    elseif hasproperty(state, :H_x) && !all(isfinite, state.H_x)
         TerminationCode.HessianNotFinite
     else
         TerminationCode.NotImplemented

src/multivariate/solvers/constrained/fminbox.jl

@@ -1,5 +1,5 @@
 using NLSolversBase:
-    value, value!, value!!, gradient, gradient!, value_gradient!, value_gradient!!
+    value, value!, gradient, gradient!, value_gradient!
 ####### FIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIX THE MIDDLE OF BOX CASE THAT WAS THERE
 mutable struct BarrierWrapper{TO,TB,Tm,TF,TDF} <: AbstractObjective
     obj::TO
@@ -69,29 +69,6 @@ function _barrier_term_gradient(x::T, l, u) where {T}
 end
 
 # Wrappers
-function NLSolversBase.value!!(bw::BarrierWrapper, x)
-    bw.Fb = _barrier_value(bw.b, x)
-    if in_box(bw, x)
-        F = value!!(bw.obj, x)
-        bw.Ftotal = muladd(bw.mu, bw.Fb, F)
-    else
-        bw.Ftotal = bw.mu * bw.Fb
-    end
-    return bw.Ftotal
-end
-function NLSolversBase.value_gradient!!(bw::BarrierWrapper, x)
-    bw.Fb = _barrier_value(bw.b, x)
-    bw.DFb .= _barrier_term_gradient.(x, bw.b.lower, bw.b.upper)
-    if in_box(bw, x)
-        F, DF = value_gradient!!(bw.obj, x)
-        bw.Ftotal = muladd(bw.mu, bw.Fb, F)
-        bw.DFtotal .= muladd.(bw.mu, bw.DFb, DF)
-    else
-        bw.Ftotal = bw.mu * bw.Fb
-        bw.DFtotal .= bw.mu .* bw.DFb
-    end
-    return bw.Ftotal, bw.DFtotal
-end
 function NLSolversBase.value_gradient!(bb::BarrierWrapper, x)
     bb.DFb .= _barrier_term_gradient.(x, bb.b.lower, bb.b.upper)
     bb.Fb = _barrier_value(bb.b, x)
@@ -115,7 +92,6 @@ function NLSolversBase.value!(obj::BarrierWrapper, x)
     end
     return obj.Ftotal
 end
-NLSolversBase.value(obj::BarrierWrapper) = obj.Ftotal
 function NLSolversBase.value(obj::BarrierWrapper, x)
     Fb = _barrier_value(obj.b, x)
     if in_box(obj, x)
@@ -134,11 +110,10 @@ function NLSolversBase.gradient!(obj::BarrierWrapper, x)
     end
     return obj.DFtotal
 end
-NLSolversBase.gradient(obj::BarrierWrapper) = obj.DFtotal
 
 # this mutates mu but not the gradients
 # Super unsafe in that it depends on x_df being correct!
-function initial_mu(obj::BarrierWrapper, F)
+function initial_mu!(obj::BarrierWrapper, g_x, F)
     T = typeof(obj.Fb) # this will not work if F is real, G is complex
     gbarrier = map(
         x ->
@@ -148,8 +123,9 @@ function initial_mu(obj::BarrierWrapper, F)
     )
 
     # obj.mu = initial_mu(gradient(obj.obj), gradient(obj.b, obj.DFb, obj.obj.x_df), T(F.mufactor), T(F.mu0))
-    obj.mu = initial_mu(gradient(obj.obj), gbarrier, T(F.mufactor), T(F.mu0))
+    return obj.mu = initial_mu(g_x, gbarrier, T(F.mufactor), T(F.mu0))
 end
+
 # Attempt to compute a reasonable default mu: at the starting
 # position, the gradient of the input function should dominate the
 # gradient of the barrier.
@@ -279,6 +255,13 @@ barrier_method(
     precondprep,
 ) = m # use `m` as is
 
+struct BoxState{T,Tx} <: ZerothOrderState
+    x::Tx
+    f_x::T
+    x_previous::Tx
+    f_x_previous::T
+end
+
 function optimize(
     f,
     l::AbstractArray,
@@ -483,14 +466,14 @@ function optimize(
     # barrier-aware optimization method instance (precondition relevance)
     _optimizer = barrier_method(F.method, P, (P, x) -> F.precondprep(P, x, l, u, dfbox))
 
+    # we wait until state has been initialized to set the initial mu because we need the gradient of the objective
     state = initial_state(_optimizer, options, dfbox, x)
-    # we wait until state has been initialized to set the initial mu because
-    # we need the gradient of the objective and initial_state will value_gradient!!
-    # the objective, so that forces an evaluation
-    if F.method isa NelderMead
+    g_x = if hasproperty(state, :g_x)
+        state.g_x
+    else
         gradient!(dfbox, x)
     end
-    dfbox.mu = initial_mu(dfbox, F)
+    initial_mu!(dfbox, g_x, F)
     if F.method isa NelderMead
         for i = 1:length(state.f_simplex)
             x = state.simplex[i]
@@ -514,10 +497,9 @@ function optimize(
     iteration = 1
 
     # define the function (dfbox) to optimize by the inner optimizer
-
     xold = copy(x)
+    fval0 = state.f_x
     _time = time()
-    fval0 = dfbox.obj.F
 
     # Optimize with current setting of mu
     if show_trace > 0
@@ -530,9 +512,7 @@ function optimize(
         println("(numbers below include barrier contribution)")
     end
 
-    # we need to update the +mu*barrier_grad part. Since we're using the
-    # value_gradient! not !! as in initial_state, we won't make a superfluous
-    # evaluation
+    # we need to update the +mu*barrier_grad part.
     if !(F.method isa NelderMead)
         value_gradient!(dfbox, x)
         reset!(_optimizer, state, dfbox, x)
@@ -681,14 +661,13 @@ function optimize(
         f_converged = _f_converged,
         g_converged = _g_converged,
     )
-    box_state = (; x, x_previous = xold, f_x_previous = fval0)
-    termination_code = _termination_code(df, g_residual(g), box_state, stopped_by, options)
+    termination_code = _termination_code(df, g_residual(g), BoxState(minimizer(results), minimum(results), xold, fval0), stopped_by, options)
 
     return MultivariateOptimizationResults(
         F,
         initial_x,
         minimizer(results),
-        df.f(minimizer(results)),
+        minimum(results),
         iteration,
         results.x_abstol,
         results.x_reltol,
@@ -699,7 +678,7 @@ function optimize(
         f_abschange(minimum(results), fval0),
         f_relchange(minimum(results), fval0),
         results.g_abstol,
-        g_residual(g, Inf),
+        g_residual(g),
         results.trace,
         results.f_calls,
         results.g_calls,