Fixed errors after merge + wrote dot() function taking a double* as second parameter

Author: cvanaret

uno/ingredients/hessian_models/HessianModelFactory.cpp

@@ -14,8 +14,7 @@
 #include "options/Options.hpp"
 
 namespace uno {
-   std::unique_ptr<HessianModel> HessianModelFactory::create([[maybe_unused]] std::optional<double> fixed_objective_multiplier,
-         const Options& options) {
+   std::unique_ptr<HessianModel> HessianModelFactory::create(double objective_multiplier, const Options& options) {
       const std::string& hessian_model = options.get_string("hessian_model");
       if (hessian_model == "exact") {
          return std::make_unique<ExactHessian>();
@@ -25,7 +24,7 @@ namespace uno {
       }
 #ifdef HAS_LAPACK
       else if (hessian_model == "L-BFGS") {
-         return std::make_unique<LBFGSHessian>(fixed_objective_multiplier, options);
+         return std::make_unique<LBFGSHessian>(objective_multiplier, options);
       }
 #endif
       else if (hessian_model == "zero") {

uno/ingredients/hessian_models/HessianModelFactory.hpp

@@ -5,7 +5,6 @@
 #define UNO_HESSIANMODELFACTORY_H
 
 #include <memory>
-#include <optional>
 
 namespace uno {
    // forward declarations
@@ -15,7 +14,7 @@ namespace uno {
 
    class HessianModelFactory {
    public:
-      static std::unique_ptr<HessianModel> create(std::optional<double> fixed_objective_multiplier, const Options& options);
+      static std::unique_ptr<HessianModel> create(double fixed_objective_multiplier, const Options& options);
    };
 } // namespace
 

uno/ingredients/hessian_models/quasi_newton/LBFGSHessian.cpp

@@ -7,18 +7,16 @@
 #include "linear_algebra/LAPACK.hpp"
 #include "linear_algebra/SymmetricMatrix.hpp"
 #include "optimization/Iterate.hpp"
+#include "optimization/OptimizationProblem.hpp"
 #include "symbolic/Expression.hpp"
 #include "symbolic/Range.hpp"
 #include "tools/Statistics.hpp"
 
 namespace uno {
-   LBFGSHessian::LBFGSHessian(std::optional<double> fixed_objective_multiplier, const Options& options):
+   LBFGSHessian::LBFGSHessian(double objective_multiplier, const Options& options):
          HessianModel(),
-         fixed_objective_multiplier(fixed_objective_multiplier),
+         objective_multiplier(objective_multiplier),
          memory_size(options.get_unsigned_int("quasi_newton_memory_size")) {
-      if (fixed_objective_multiplier.has_value()) {
-         DEBUG << "L-BFGS Hessian model was declared with a fixed objective multiplier of " << *fixed_objective_multiplier << '\n';
-      }
    }
 
    bool LBFGSHessian::has_implicit_representation() const {
@@ -67,8 +65,11 @@ namespace uno {
 
    // Hessian-vector product where the Hessian approximation is Bk = B0 - U U^T + V V^T and B0 = delta I
    // Bk v = (B0 - U U^T + V V^T) v = delta v - U U^T x + V V^T x
-   void LBFGSHessian::compute_hessian_vector_product(const Model& model, const double* vector, double /*objective_multiplier*/,
+   void LBFGSHessian::compute_hessian_vector_product(const Model& model, const double* vector, double objective_multiplier,
          const Vector<double>& /*constraint_multipliers*/, double* result) {
+      assert(objective_multiplier == this->objective_multiplier &&
+         "The L-BFGS Hessian model was initialized with a different objective multiplier");
+
       if (this->hessian_recomputation_required) {
          this->recompute_hessian_representation();
          this->hessian_recomputation_required = false;
@@ -136,23 +137,21 @@ namespace uno {
    // fill the Y matrix: y = \nabla L(x_k, y_k, z_k) - \nabla L(x_{k-1}, y_k, z_k)
    void LBFGSHessian::update_Y_matrix(const Model& model, Iterate& current_iterate, Iterate& trial_iterate) {
       // evaluate Lagrangian gradients at the current and trial iterates, both with the trial multipliers
-      // TODO objective multiplier is hardcoded for the moment
-      if (this->fixed_objective_multiplier.has_value()) {
-         current_iterate.evaluate_objective_gradient(model);
-         current_iterate.evaluate_constraint_jacobian(model);
-         trial_iterate.evaluate_objective_gradient(model);
-         trial_iterate.evaluate_constraint_jacobian(model);
-         // TODO preallocate
-         Vector<double> current_lagrangian_gradient(this->dimension);
-         Vector<double> trial_lagrangian_gradient(this->dimension);
-         const double objective_multiplier = *this->fixed_objective_multiplier;
-         model.evaluate_lagrangian_gradient(current_lagrangian_gradient, current_iterate, trial_iterate.multipliers, objective_multiplier);
-         model.evaluate_lagrangian_gradient(trial_lagrangian_gradient, trial_iterate, trial_iterate.multipliers, objective_multiplier);
-         this->Y_matrix.column(this->current_memory_slot) = trial_lagrangian_gradient - current_lagrangian_gradient;
-      }
-      else {
-         throw std::runtime_error("LBFGSHessian::update_Y_matrix: the objective multiplier varies. This is not implemented yet");
-      }
+      current_iterate.evaluate_objective_gradient(model);
+      current_iterate.evaluate_constraint_jacobian(model);
+      trial_iterate.evaluate_objective_gradient(model);
+      trial_iterate.evaluate_constraint_jacobian(model);
+      // TODO preallocate
+      LagrangianGradient<double> current_split_lagrangian_gradient(this->dimension);
+      LagrangianGradient<double> trial_split_lagrangian_gradient(this->dimension);
+      const OptimizationProblem problem{model};
+      problem.evaluate_lagrangian_gradient(current_split_lagrangian_gradient, current_iterate, trial_iterate.multipliers);
+      problem.evaluate_lagrangian_gradient(trial_split_lagrangian_gradient, trial_iterate, trial_iterate.multipliers);
+      const auto current_lagrangian_gradient = this->objective_multiplier * current_split_lagrangian_gradient.objective_contribution
+         + current_split_lagrangian_gradient.constraints_contribution;
+      const auto trial_lagrangian_gradient = this->objective_multiplier * trial_split_lagrangian_gradient.objective_contribution
+         + trial_split_lagrangian_gradient.constraints_contribution;
+      this->Y_matrix.column(this->current_memory_slot) = trial_lagrangian_gradient - current_lagrangian_gradient;
    }
 
    void LBFGSHessian::update_D_matrix() {

uno/ingredients/hessian_models/quasi_newton/LBFGSHessian.hpp

@@ -4,7 +4,6 @@
 #ifndef UNO_LBFGSHESSIAN_H
 #define UNO_LBFGSHESSIAN_H
 
-#include <optional>
 #include <vector>
 #include "../HessianModel.hpp"
 #include "linear_algebra/DenseMatrix.hpp"
@@ -19,7 +18,7 @@ namespace uno {
    // J J^T = M = Sk^T B0 Sk + Lk Dk^(-1) Lk^T
    class LBFGSHessian: public HessianModel {
    public:
-      LBFGSHessian(std::optional<double> fixed_objective_multiplier, const Options& options);
+      LBFGSHessian(double objective_multiplier, const Options& options);
       ~LBFGSHessian() override = default;
 
       [[nodiscard]] bool has_implicit_representation() const override;
@@ -38,7 +37,7 @@ namespace uno {
 
    protected:
       size_t dimension{};
-      std::optional<double> fixed_objective_multiplier;
+      const double objective_multiplier;
       const size_t memory_size; // user defined
       size_t number_entries_in_memory{0}; // 0 <= used_memory_size <= memory_size
       size_t current_memory_slot{0}; // 0 <= current_available_slot < memory_size

uno/model/Model.cpp

@@ -5,19 +5,16 @@
 #include <utility>
 #include "Model.hpp"
 #include "linear_algebra/Vector.hpp"
-#include "optimization/Iterate.hpp"
-#include "optimization/LagrangianGradient.hpp"
-#include "optimization/Multipliers.hpp"
 
 namespace uno {
+   // abstract Problem class
    Model::Model(std::string name, size_t number_variables, size_t number_constraints, double objective_sign) :
          name(std::move(name)), number_variables(number_variables), number_constraints(number_constraints), objective_sign(objective_sign) {
    }
 
    void Model::project_onto_variable_bounds(Vector<double>& x) const {
       for (size_t variable_index: Range(this->number_variables)) {
-         x[variable_index] = std::max(std::min(x[variable_index], this->variable_upper_bound(variable_index)),
-            this->variable_lower_bound(variable_index));
+         x[variable_index] = std::max(std::min(x[variable_index], this->variable_upper_bound(variable_index)), this->variable_lower_bound(variable_index));
       }
    }
 
@@ -31,54 +28,4 @@ namespace uno {
       const double upper_bound_violation = std::max(0., constraint_value - this->constraint_upper_bound(constraint_index));
       return std::max(lower_bound_violation, upper_bound_violation);
    }
-
-   void Model::evaluate_lagrangian_gradient(Vector<double>& lagrangian_gradient, const Iterate& iterate,
-         const Multipliers& multipliers, double objective_multiplier) const {
-      lagrangian_gradient.fill(0.);
-      // compute and scale the objective contribution
-      if (objective_multiplier != 0.) {
-         this->evaluate_lagrangian_gradient_objective(lagrangian_gradient, iterate);
-         lagrangian_gradient.scale(objective_multiplier);
-      }
-      // add the constraints contribution
-      this->evaluate_lagrangian_gradient_constraints(lagrangian_gradient, iterate, multipliers);
-   }
-
-   // Lagrangian gradient split in two parts: objective contribution and constraints' contribution
-   void Model::evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient, const Iterate& iterate,
-         const Multipliers& multipliers) const {
-      lagrangian_gradient.objective_contribution.fill(0.);
-      lagrangian_gradient.constraints_contribution.fill(0.);
-      this->evaluate_lagrangian_gradient_objective(lagrangian_gradient.objective_contribution, iterate);
-      this->evaluate_lagrangian_gradient_constraints(lagrangian_gradient.constraints_contribution, iterate, multipliers);
-   }
-
-   // protected member functions
-
-   // Lagrangian gradient split in two parts: objective contribution and constraints' contribution
-   void Model::evaluate_lagrangian_gradient_objective(Vector<double>& objective_contribution, const Iterate& iterate) const {
-      // objective gradient
-      /*
-      for (auto [variable_index, derivative]: iterate.evaluations.objective_gradient) {
-         objective_contribution[variable_index] += derivative;
-      }
-      */
-   }
-
-   // Lagrangian gradient split in two parts: objective contribution and constraints' contribution
-   void Model::evaluate_lagrangian_gradient_constraints(Vector<double>& constraints_contribution, const Iterate& iterate,
-         const Multipliers& multipliers) const {
-      // constraints
-      for (size_t constraint_index: Range(this->number_constraints)) {
-         if (multipliers.constraints[constraint_index] != 0.) {
-            for (auto [variable_index, derivative]: iterate.evaluations.constraint_jacobian[constraint_index]) {
-               constraints_contribution[variable_index] -= multipliers.constraints[constraint_index] * derivative;
-            }
-         }
-      }
-      // bound constraints
-      for (size_t variable_index: Range(this->number_variables)) {
-         constraints_contribution[variable_index] -= (multipliers.lower_bounds[variable_index] + multipliers.upper_bounds[variable_index]);
-      }
-   }
 } // namespace

uno/model/Model.hpp

@@ -15,9 +15,6 @@ namespace uno {
    template <typename ElementType>
    class Collection;
    template <typename ElementType>
-   class LagrangianGradient;
-   class Multipliers;
-   template <typename ElementType>
    class RectangularMatrix;
    template <typename ElementType>
    class SparseVector;
@@ -88,16 +85,6 @@ namespace uno {
       [[nodiscard]] virtual double constraint_violation(double constraint_value, size_t constraint_index) const;
       template <typename Array>
       double constraint_violation(const Array& constraints, Norm residual_norm) const;
-
-      void evaluate_lagrangian_gradient(Vector<double>& lagrangian_gradient, const Iterate& iterate,
-            const Multipliers& multipliers, double objective_multiplier) const;
-      void evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient, const Iterate& iterate,
-         const Multipliers& multipliers) const;
-
-   protected:
-      void evaluate_lagrangian_gradient_objective(Vector<double>& objective_contribution, const Iterate& iterate) const;
-      void evaluate_lagrangian_gradient_constraints(Vector<double>& constraints_contribution, const Iterate& iterate,
-         const Multipliers& multipliers) const;
    };
 
    // compute ||c||
@@ -111,4 +98,4 @@ namespace uno {
    }
 } // namespace
 
-#endif // UNO_MODEL_H
+#endif // UNO_MODEL_H

uno/optimization/LagrangianGradient.hpp

@@ -8,57 +8,6 @@
 #include "linear_algebra/Vector.hpp"
 
 namespace uno {
-   // forward declaration
-   template <typename ElementType>
-   class LagrangianGradient;
-
-   // The LagrangianGradient class represents the dense Lagrangian gradient broken into:
-   // - the objective contribution \nabla f(x_k)
-   // - the constraint contribution -\nabla c(x_k) y_k - z_k
-   // The two contributions can be assembled with a given objective multiplier \rho (see assemble(double)):
-   // \rho \nabla f(x_k) - \nabla c(x_k) y_k - z_k
-   // The resulting object, an AssembledLagrangianGradient, is a wrapper around the LagrangianGradient and \rho
-
-   template <typename ElementType>
-   class AssembledLagrangianGradient {
-   public:
-      using value_type = ElementType;
-
-      AssembledLagrangianGradient(const LagrangianGradient<ElementType>& lagrangian_gradient, ElementType objective_multiplier);
-
-      [[nodiscard]] size_t size() const;
-      [[nodiscard]] ElementType operator[](size_t variable_index) const;
-
-   protected:
-      const LagrangianGradient<ElementType>& lagrangian_gradient;
-      const ElementType objective_multiplier;
-   };
-
-   template<typename ElementType>
-   AssembledLagrangianGradient<ElementType>::AssembledLagrangianGradient(const LagrangianGradient<ElementType> &lagrangian_gradient,
-      ElementType objective_multiplier): lagrangian_gradient(lagrangian_gradient), objective_multiplier(objective_multiplier) { }
-
-   template <typename ElementType>
-   size_t AssembledLagrangianGradient<ElementType>::size() const {
-      return this->lagrangian_gradient.size();
-   }
-
-  // access i-th element
-  template <typename ElementType>
-  ElementType AssembledLagrangianGradient<ElementType>::operator[](size_t variable_index) const {
-     return this->objective_multiplier * this->lagrangian_gradient.objective_contribution[variable_index] +
-        this->lagrangian_gradient.constraints_contribution[variable_index];
-  }
-
-   template <typename ElementType>
-   std::ostream& operator<<(std::ostream& stream, const AssembledLagrangianGradient<ElementType>& gradient) {
-      for (size_t variable_index: Range(gradient.size())) {
-         stream << gradient[variable_index] << ' ';
-      }
-      stream << '\n';
-      return stream;
-   }
-
    // Gradient of the Lagrangian
    // Keep the objective and constraint contributions separate. This helps:
    // - computing the KKT and FJ stationarity conditions
@@ -73,8 +22,8 @@ namespace uno {
 
       explicit LagrangianGradient(size_t number_variables);
       [[nodiscard]] size_t size() const;
+      [[nodiscard]] ElementType operator[](size_t variable_index) const;
       void resize(size_t number_variables);
-      AssembledLagrangianGradient<ElementType> assemble(double objective_multiplier) const;
    };
 
    template <typename ElementType>
@@ -88,16 +37,25 @@ namespace uno {
       return this->objective_contribution.size();
    }
 
+   // access i-th element
+   template <typename ElementType>
+   ElementType LagrangianGradient<ElementType>::operator[](size_t variable_index) const {
+      return this->objective_contribution[variable_index] + this->constraints_contribution[variable_index];
+   }
+
    template <typename ElementType>
    void LagrangianGradient<ElementType>::resize(size_t number_variables) {
       this->objective_contribution.resize(number_variables);
       this->constraints_contribution.resize(number_variables);
    }
 
    template <typename ElementType>
-   AssembledLagrangianGradient<ElementType> LagrangianGradient<ElementType>::assemble(double objective_multiplier) const {
-      // TODO use existing Expression
-      return {*this, objective_multiplier};
+   std::ostream& operator<<(std::ostream& stream, const LagrangianGradient<ElementType>& gradient) {
+      for (size_t variable_index: Range(gradient.constraints_contribution.size())) {
+         stream << gradient[variable_index] << ' ';
+      }
+      stream << '\n';
+      return stream;
    }
 } // namespace