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

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,28 +28,19 @@ 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)
-    gradient!(obj.inner_obj, xin)
-    project_tangent!(obj.manifold, gradient(obj.inner_obj), xin)
-    return gradient(obj.inner_obj)
+    g_xin = gradient!(obj.inner_obj, xin)
+    project_tangent!(obj.manifold, g_xin, xin)
+    return g_xin
 end
 function NLSolversBase.value_gradient!(obj::ManifoldObjective, x)
     xin = retract(obj.manifold, x)
-    value_gradient!(obj.inner_obj, xin)
-    project_tangent!(obj.manifold, gradient(obj.inner_obj), xin)
-    return value(obj.inner_obj)
+    f_xin, g_xin = value_gradient!(obj.inner_obj, xin)
+    project_tangent!(obj.manifold, g_xin, xin)
+    return f_xin, g_xin
 end
 
 """Flat Euclidean space {R,C}^N, with projections equal to the identity."""

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,21 +1,37 @@
-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
@@ -24,14 +40,14 @@ update_h!(d, state, method::SecondOrderOptimizer) = hessian!(d, state.x)
 
 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 +57,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 +67,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 +129,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 +143,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 +157,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 +170,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 +202,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 +226,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

@@ -72,64 +72,75 @@ function gradient(bb::BoxBarrier, g, x)
     g = copy(g)
     g .= _barrier_term_gradient.(x, bb.lower, bb.upper)
 end
+
 # Wrappers
 function value!!(bw::BarrierWrapper, x)
     bw.Fb = value(bw.b, x)
-    bw.Ftotal = bw.mu * bw.Fb
     if in_box(bw, x)
-        value!!(bw.obj, x)
-        bw.Ftotal += value(bw.obj)
+        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 value_gradient!!(bw::BarrierWrapper, x)
     bw.Fb = value(bw.b, x)
-    bw.Ftotal = bw.mu * bw.Fb
     bw.DFb .= _barrier_term_gradient.(x, bw.b.lower, bw.b.upper)
-    bw.DFtotal .= bw.mu .* bw.DFb
-    if in_box(bw, x)
-        value_gradient!!(bw.obj, x)
-        bw.Ftotal += value(bw.obj)
-        bw.DFtotal .+= gradient(bw.obj)
+    if in_box(bw.b, 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 value_gradient!(bb::BarrierWrapper, x)
     bb.DFb .= _barrier_term_gradient.(x, bb.b.lower, bb.b.upper)
     bb.Fb = value(bb.b, x)
-    bb.DFtotal .= bb.mu .* bb.DFb
-    bb.Ftotal = bb.mu * bb.Fb
-
-    if in_box(bb, x)
-        value_gradient!(bb.obj, x)
-        bb.DFtotal .+= gradient(bb.obj)
-        bb.Ftotal += value(bb.obj)
+    if in_box(bb.b, x)
+        F, DF = value_gradient!(bb.obj, x)
+        bb.DFtotal .= muladd.(bb.mu, bb.DFb, DF)
+        bb.Ftotal = muladd(bb.mu, bb.Fb, F)
+    else
+        bb.DFtotal .= bb.mu .* bb.DFb
+        bb.Ftotal = bb.mu * bb.Fb
     end
+    return bb.Ftotal, bb.DFtotal
 end
 value(bb::BoxBarrier, x) =
     mapreduce(x -> _barrier_term_value(x...), +, zip(x, bb.lower, bb.upper))
 function value!(obj::BarrierWrapper, x)
     obj.Fb = value(obj.b, x)
-    obj.Ftotal = obj.mu * obj.Fb
-    if in_box(obj, x)
-        value!(obj.obj, x)
-        obj.Ftotal += value(obj.obj)
+    if in_box(obj.b, x)
+        F = value!(obj.obj, x)
+        obj.Ftotal = muladd(obj.mu, obj.Fb, F)
+    else
+        obj.Ftotal = obj.mu * obj.Fb
     end
-    obj.Ftotal
+    return obj.Ftotal
 end
-value(obj::BarrierWrapper) = obj.Ftotal
+
 function value(obj::BarrierWrapper, x)
-    F = obj.mu * value(obj.b, x)
-    if in_box(obj, x)
-        F += value(obj.obj, x)
+    Fb = value(obj.b, x)
+    if in_box(obj.b, x)
+        return muladd(obj.mu, Fb, value(obj.obj, x))
+    else
+        return obj.mu * Fb
     end
-    F
 end
 function gradient!(obj::BarrierWrapper, x)
-    gradient!(obj.obj, x)
-    obj.DFb .= gradient(obj.b, obj.DFb, x) # this should just be inplace?
-    obj.DFtotal .= gradient(obj.obj) .+ obj.mu * obj.Fb
+    obj.DFb .= _barrier_term_gradient.(x, obj.b.lower, obj.b.upper)
+    if in_box(obj.b, x)
+        DF = gradient!(obj.obj, x)
+        obj.DFtotal .= muladd.(obj.mu, obj.Fb, DF)
+    else
+        obj.DFtotal .= obj.mu .* obj.DFb
+    end
+    return obj.DFtotal
 end
-gradient(obj::BarrierWrapper) = obj.DFtotal
 
 # this mutates mu but not the gradients
 # Super unsafe in that it depends on x_df being correct!
@@ -694,7 +705,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,