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

Author: devmotion

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/multivariate/optimize/optimize.jl

@@ -1,21 +1,38 @@
 update_g!(d, state, method) = 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)
     # 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::ZerothOrderOptimizer)
+    f_x = value!(d, state.x)
+    state.f_x = f_x
+    return nothing
+end
 function update_fg!(d, state, method::FirstOrderOptimizer)
-    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_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 +41,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,
@@ -51,7 +68,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 +130,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
@@ -141,7 +158,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),
@@ -162,10 +179,10 @@ function optimize(
         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 +203,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 +227,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

@@ -77,59 +77,66 @@ 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 += F
     end
 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)
+        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)
+        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)
+        F = value!(obj.obj, x)
+        obj.Ftotal += F
     end
     obj.Ftotal
 end
-value(obj::BarrierWrapper) = obj.Ftotal
+
 function value(obj::BarrierWrapper, x)
-    F = obj.mu * value(obj.b, x)
+    Fb = value(obj.b, x)
     if in_box(obj, x)
-        F += value(obj.obj, 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!

src/multivariate/solvers/constrained/ipnewton/interior.jl

@@ -209,9 +209,9 @@ ls_update!(
 function initial_convergence(d, state, method::ConstrainedOptimizer, initial_x, options)
     # TODO: Make sure state.bgrad has been evaluated at initial_x
     # state.bgrad normally comes from constraints.c!(..., initial_x) in initial_state
-    gradient!(d, initial_x)
-    stopped = !isfinite(value(d)) || any(!isfinite, gradient(d))
-    g_residual(d, state) + norm(state.bgrad, Inf) < options.g_abstol, stopped
+    f_x, g_x = gradient!(d, initial_x)
+    stopped = !isfinite(f_x) || any(!isfinite, g_x)
+    g_residual(g_x, state) + norm(state.bgrad, Inf) < options.g_abstol, stopped
 end
 
 function optimize(
@@ -342,10 +342,10 @@ function optimize(
         x_relchange(state),
         T(options.f_abstol),
         T(options.f_reltol),
-        f_abschange(d, state),
-        f_relchange(d, state),
+        f_abschange(state),
+        f_relchange(state),
         T(options.g_abstol),
-        g_residual(d, state),
+        g_residual(state),
         tr,
         f_calls(d),
         g_calls(d),

src/multivariate/solvers/constrained/ipnewton/utilities/trace.jl

@@ -45,7 +45,7 @@ function trace!(tr, d, state, iteration, method::IPOptimizer, options, curr_time
     update!(
         tr,
         iteration,
-        value(d),
+        state.f_x,
         g_norm,
         dt,
         options.store_trace,

src/multivariate/solvers/constrained/samin.jl

@@ -225,10 +225,10 @@ function optimize(
                             NaN,# x_abschange(state),
                             f_tol,#T(options.f_tol),
                             0.0,#T(options.f_tol),
-                            f_absΔ,#f_abschange(d, state),
-                            NaN,#f_abschange(d, state),
+                            f_absΔ,#f_abschange(state),
+                            NaN,#f_abschange(state),
                             0.0,#T(options.g_tol),
-                            NaN,#g_residual(d),
+                            NaN,#g_residual(state),
                             tr,
                             f_calls(d),
                             g_calls(d),

src/multivariate/solvers/first_order/accelerated_gradient_descent.jl

@@ -35,21 +35,21 @@ end
 
 function initial_state(
     method::AcceleratedGradientDescent,
-    options,
+    options::Options,
     d,
-    initial_x::AbstractArray{T},
-) where {T}
+    initial_x::AbstractArray,
+)
     initial_x = copy(initial_x)
     retract!(method.manifold, initial_x)
-
-    value_gradient!!(d, initial_x)
-
-    project_tangent!(method.manifold, gradient(d), initial_x)
+    f_x, g_x = value_gradient!(d, initial_x)
+    project_tangent!(method.manifold, g_x, initial_x)
 
     AcceleratedGradientDescentState(
         copy(initial_x), # Maintain current state in state.x
-        copy(initial_x), # Maintain previous state in state.x_previous
-        real(T)(NaN), # Store previous f in state.f_x_previous
+        copy(g_x), # Maintain current gradient in state.g_x
+        f_x, # Maintain current f in state.f_x
+        fill!(similar(initial_x), NaN), # Maintain previous state in state.x_previous
+        oftype(f_x, NaN), # Store previous f in state.f_x_previous
         0, # Iteration
         copy(initial_x), # Maintain intermediary current state in state.y
         similar(initial_x), # Maintain intermediary state in state.y_previous
@@ -63,11 +63,14 @@ function update_state!(
     state::AcceleratedGradientDescentState,
     method::AcceleratedGradientDescent,
 )
-    value_gradient!(d, state.x)
+    f_x, g_x = value_gradient!(d, state.x)
     state.iteration += 1
-    project_tangent!(method.manifold, gradient(d), state.x)
+    project_tangent!(method.manifold, g_x, state.x)
+    copyto!(state.g_x, g_x)
+    state.f_x = f_x
+
     # Search direction is always the negative gradient
-    state.s .= .-gradient(d)
+    state.s .= .-g_x
 
     # Determine the distance of movement along the search line
     lssuccess = perform_linesearch!(state, method, ManifoldObjective(method.manifold, d))