Reduce amount of work in State

Author: PatWie

.github/workflows/ci.yml

@@ -58,3 +58,4 @@ jobs:
       - name: Build
         run: |
           bazel run simple
+          bazel run constrained_simple

README.md

@@ -103,18 +103,21 @@ int main(int argc, char const *argv[]) {
     std::cout << "Initial point: " << x.transpose() << std::endl;
 
     // Evaluate
-    auto state = f.GetState(x);
-    std::cout << "Function value at initial point: " << f(x) << std::endl;
-    std::cout << "Gradient at initial point: " << state.gradient << std::endl;
+    Eigen::VectorXd gradient(2);
+    double value = f(x, &gradient);
+    std::cout << "Function value at initial point: " << value << std::endl;
+    std::cout << "Gradient at initial point: " << gradient << std::endl;
 
     // Minimize the Rosenbrock function using the BFGS solver.
     using Solver = cppoptlib::solver::Bfgs<Rosenbrock>;
+    auto init_state = f.GetState(x);
     Solver solver;
-    auto [solution_state, solver_progress] = solver.Minimize(f, x);
+    auto [solution_state, solver_progress] = solver.Minimize(f, init_state);
 
     // Display the results of the optimization.
     std::cout << "Optimal solution: " << solution_state.x.transpose() << std::endl;
-    std::cout << "Optimal function value: " << solution_state.value << std::endl;
+    std::cout << "Optimal function value: " << f(solution_state.x) << std::endl;
+
     std::cout << "Number of iterations: " << solver_progress.num_iterations << std::endl;
     std::cout << "Solver status: " << solver_progress.status << std::endl;
 

include/cppoptlib/constrained_function.h

@@ -23,10 +23,7 @@ struct ConstrainedState : public State<base_t> {
   using unconstrained_base_t = FunctionBase<typename base_t::scalar_t,
                                             base_t::Dim, base_t::DiffLevel, 0>;
 
-  typename base_t::scalar_t value = 0;
   typename base_t::vector_t x;
-  typename base_t::vector_t gradient;
-
   std::array<typename base_t::scalar_t, base_t::NumConstraints>
       lagrange_multipliers;
   std::array<typename base_t::scalar_t, base_t::NumConstraints> violations;
@@ -46,19 +43,15 @@ struct ConstrainedState : public State<base_t> {
   }
 
   void CopyState(const self_t &rhs) {
-    value = rhs.value;
     x = rhs.x.eval();
-    gradient = rhs.gradient.eval();
     penalty = rhs.penalty;
     lagrange_multipliers = rhs.lagrange_multipliers;
     violations = rhs.violations;
   }
 
   State<unconstrained_base_t> AsUnconstrained() const {
     State<unconstrained_base_t> state;
-    state.value = value;
     state.x = x.eval();
-    state.gradient = gradient.eval();
     return state;
   }
 };
@@ -90,18 +83,7 @@ class UnconstrainedFunctionAdapter
     const typename cfunction_t::state_t inner = constrained_function.GetState(
         x, constrained_state.lagrange_multipliers, constrained_state.penalty);
     typename cfunction_t::unconstrained_function_t::state_t unconstrained_state;
-    unconstrained_state.value = inner.value;
     unconstrained_state.x = inner.x;
-    if constexpr ((cfunction_t::unconstrained_function_t::Differentiability ==
-                   Differentiability::First) ||
-                  (cfunction_t::unconstrained_function_t::Differentiability ==
-                   Differentiability::Second)) {
-      unconstrained_state.gradient = inner.gradient;
-    }
-    if constexpr (cfunction_t::unconstrained_function_t::Differentiability ==
-                  Differentiability::Second) {
-      unconstrained_state.hessian = inner.hessian;
-    }
     return unconstrained_state;
   }
 
@@ -185,30 +167,12 @@ struct ConstrainedFunction {
 
     state_t constrained_state;
     constrained_state.x = objective_state.x;
-    constrained_state.value = objective_state.value;
-    constrained_state.gradient = objective_state.gradient;
+    constrained_state.penalty = penalty;
 
-    // Sum augmented penalties for hard constraints.
     for (std::size_t i = 0; i < TNumConstraints; ++i) {
-      const typename function_t::state_t constraint_state =
-          constraints_[i]->GetState(x);
-      const scalar_t cost = constraint_state.value;
-      const scalar_t violation = cost;
-
-      const scalar_t lambda = lagrange_multipliers[i];
-      const scalar_t aug_cost =
-          violation + lambda * violation +
-          static_cast<scalar_t>(0.5) * penalty * violation * violation;
-      constrained_state.value += aug_cost;
-      // Augmented gradient (only active if the constraint is violated).
-      const scalar_t a = scalar_t(1) + lambda + penalty * violation;
-      const typename base_t::vector_t scaled_local_grad =
-          a * constraint_state.gradient;
-      typename base_t::vector_t aug_grad =
-          (cost > scalar_t(0)) ? scaled_local_grad
-                               : base_t::vector_t::Zero(x.size());
-      constrained_state.gradient += aug_grad;
+      const scalar_t violation = constraints_[i]->operator()(x);
       constrained_state.violations[i] = violation;
+      constrained_state.lagrange_multipliers[i] = lagrange_multipliers[i];
     }
 
     return constrained_state;

include/cppoptlib/function.h

@@ -45,75 +45,20 @@ class Function : public FunctionBase<TScalar, TDim, TDifferentiability> {
   using state_t = State<Function<TScalar, TDim, TDifferentiability>>;
 };
 
-template <class TScalar, int TDim>
-struct State<FunctionBase<TScalar, TDim, Differentiability::None>> {
-  using base_t = FunctionBase<TScalar, TDim, Differentiability::None>;
-  using state_t = State<base_t>;
-
-  typename base_t::scalar_t value = 0;
-  typename base_t::vector_t x;
-
-  State() = default;
-
-  State(const state_t &rhs) : value(rhs.value), x(rhs.x.eval()) {}
-
-  state_t &operator=(const state_t &rhs) {
-    if (this != &rhs) {
-      value = rhs.value;
-      x = rhs.x.eval();
-    }
-    return *this;
-  }
-};
-
-template <class TScalar, int TDim>
-struct State<FunctionBase<TScalar, TDim, Differentiability::First>> {
-  using base_t = FunctionBase<TScalar, TDim, Differentiability::First>;
-  using state_t = State<base_t>;
-
-  typename base_t::scalar_t value = 0;
-  typename base_t::vector_t x;
-  typename base_t::vector_t gradient;
-
-  State() = default;
-
-  State(const state_t &rhs)
-      : value(rhs.value), x(rhs.x.eval()), gradient(rhs.gradient.eval()) {}
-
-  state_t &operator=(const state_t &rhs) {
-    if (this != &rhs) {
-      value = rhs.value;
-      x = rhs.x.eval();
-      gradient = rhs.gradient.eval();
-    }
-    return *this;
-  }
-};
-
-template <class TScalar, int TDim>
-struct State<FunctionBase<TScalar, TDim, Differentiability::Second>> {
-  using base_t = FunctionBase<TScalar, TDim, Differentiability::Second>;
+template <class TScalar, int TDim, Differentiability TDifferentiability>
+struct State<FunctionBase<TScalar, TDim, TDifferentiability>> {
+  using base_t = FunctionBase<TScalar, TDim, TDifferentiability>;
   using state_t = State<base_t>;
 
-  typename base_t::scalar_t value = 0;
   typename base_t::vector_t x;
-  typename base_t::vector_t gradient;
-  typename base_t::matrix_t hessian;
 
   State() = default;
 
-  State(const state_t &rhs)
-      : value(rhs.value),
-        x(rhs.x.eval()),
-        gradient(rhs.gradient.eval()),
-        hessian(rhs.hessian.eval()) {}
+  State(const state_t &rhs) : x(rhs.x.eval()) {}
 
   state_t &operator=(const state_t &rhs) {
     if (this != &rhs) {
-      value = rhs.value;
       x = rhs.x.eval();
-      gradient = rhs.gradient.eval();
-      hessian = rhs.hessian.eval();
     }
     return *this;
   }
@@ -141,7 +86,6 @@ class Function<TScalar, TDim, Differentiability::None>
   state_t GetState(const typename base_t::vector_t &x) const {
     state_t state;
     state.x = x;
-    state.value = this->operator()(x);
     return state;
   }
 };
@@ -168,7 +112,6 @@ class Function<TScalar, TDim, Differentiability::First>
   state_t GetState(const typename base_t::vector_t &x) const {
     state_t state;
     state.x = x;
-    state.value = this->operator()(x, &state.gradient);
     return state;
   }
 };
@@ -196,7 +139,6 @@ class Function<TScalar, TDim, Differentiability::Second>
   state_t GetState(const typename base_t::vector_t &x) const {
     state_t state;
     state.x = x;
-    state.value = this->operator()(x, &state.gradient, &state.hessian);
     return state;
   }
 };

include/cppoptlib/solver/augmented_lagrangian.h

@@ -34,7 +34,8 @@ class AugmentedLagrangian : public Solver<function_t> {
   AugmentedLagrangian(const solver_t &inner_solver)
       : inner_solver_(inner_solver) {}
 
-  void InitializeSolver(const state_t & /*initial_state*/) override {}
+  void InitializeSolver(const function_t & /*function*/,
+                        const state_t & /*initial_state*/) override {}
 
   state_t OptimizationStep(const function_t &function, const state_t &state,
                            const progress_t & /*progress*/) override {

include/cppoptlib/solver/bfgs.h

@@ -32,32 +32,38 @@ class Bfgs : public Solver<function_t> {
 
   using Superclass::Superclass;
 
-  void InitializeSolver(const state_t &initial_state) override {
+  void InitializeSolver(const function_t & /*function*/,
+                        const state_t &initial_state) override {
     dim_ = initial_state.x.rows();
     inverse_hessian_ =
         matrix_t::Identity(initial_state.x.rows(), initial_state.x.rows());
   }
 
   state_t OptimizationStep(const function_t &function, const state_t &current,
                            const progress_t & /*progress*/) override {
-    vector_t search_direction = -inverse_hessian_ * current.gradient;
+    vector_t current_gradient;
+    function(current.x, &current_gradient);
+
+    vector_t search_direction = -inverse_hessian_ * current_gradient;
 
     // If not positive definit re-initialize Hessian.
-    const scalar_t phi = current.gradient.dot(search_direction);
+    const scalar_t phi = current_gradient.dot(search_direction);
     if ((phi > 0) || std::isnan(phi)) {
       // no, we reset the hessian approximation
       inverse_hessian_ = matrix_t::Identity(dim_, dim_);
-      search_direction = -current.gradient;
+      search_direction = -current_gradient;
     }
 
     const scalar_t rate = linesearch::MoreThuente<function_t, 1>::Search(
         current.x, search_direction, function);
 
     const state_t next = function.GetState(current.x + rate * search_direction);
+    vector_t next_gradient;
+    function(next.x, &next_gradient);
 
     // Update inverse Hessian estimate.
     const vector_t s = rate * search_direction;
-    const vector_t y = next.gradient - current.gradient;
+    const vector_t y = next_gradient - current_gradient;
     const scalar_t rho = 1.0 / y.dot(s);
 
     inverse_hessian_ =

include/cppoptlib/solver/conjugated_gradient_descent.h

@@ -32,20 +32,23 @@ class ConjugatedGradientDescent : public Solver<function_t> {
 
   using Superclass::Superclass;
 
-  void InitializeSolver(const state_t &initial_state) override {
-    previous_gradient_ = initial_state.gradient;
+  void InitializeSolver(const function_t &function,
+                        const state_t &initial_state) override {
+    function(initial_state.x, &previous_gradient_);
   }
 
   state_t OptimizationStep(const function_t &function, const state_t &current,
                            const progress_t &progress) override {
+    vector_t current_gradient;
+    function(current.x, &current_gradient);
     if (progress.num_iterations == 0) {
-      search_direction_ = -current.gradient;
+      search_direction_ = -current_gradient;
     } else {
-      const double beta = current.gradient.dot(current.gradient) /
+      const double beta = current_gradient.dot(current_gradient) /
                           (previous_gradient_.dot(previous_gradient_));
-      search_direction_ = -current.gradient + beta * search_direction_;
+      search_direction_ = -current_gradient + beta * search_direction_;
     }
-    previous_gradient_ = current.gradient;
+    previous_gradient_ = current_gradient;
 
     const scalar_t rate = linesearch::Armijo<function_t, 1>::Search(
         current.x, search_direction_, function);

include/cppoptlib/solver/gradient_descent.h

@@ -31,14 +31,17 @@ class GradientDescent : public Solver<function_t> {
 
   using Superclass::Superclass;
 
-  void InitializeSolver(const state_t & /*initial_state*/) override {}
+  void InitializeSolver(const function_t & /*function*/,
+                        const state_t & /*initial_state*/) override {}
 
   state_t OptimizationStep(const function_t &function, const state_t &current,
                            const progress_t & /*progress*/) override {
+    vector_t gradient;
+    function(current.x, &gradient);
     const scalar_t rate = linesearch::MoreThuente<function_t, 1>::Search(
-        current.x, -current.gradient, function);
+        current.x, -gradient, function);
 
-    return function.GetState(current.x - rate * current.gradient);
+    return function.GetState(current.x - rate * gradient);
   }
 };