Merge pull request #128 from cvanaret/refactoring_qp

Allocation and creation of QP/LP within QP/LP solver

uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp

@@ -72,24 +72,15 @@ namespace uno {
    }
 
    std::function<double(double)> ConstraintRelaxationStrategy::compute_predicted_objective_reduction_model(const Iterate& current_iterate,
-         const Vector<double>& primal_direction, double step_length, const SymmetricMatrix<size_t, double>& hessian) const {
+         const Vector<double>& primal_direction, double step_length) const {
       // predicted objective reduction: "-∇f(x)^T (αd) - α^2/2 d^T H d"
       const double directional_derivative = dot(primal_direction, current_iterate.evaluations.objective_gradient);
-      const double quadratic_term = hessian.quadratic_product(primal_direction, primal_direction);
+      const double quadratic_term = this->first_order_predicted_reduction ? 0. : this->inequality_handling_method->hessian_quadratic_product(primal_direction);
       return [=](double objective_multiplier) {
          return step_length * (-objective_multiplier*directional_derivative) - step_length*step_length/2. * quadratic_term;
       };
    }
 
-   std::function<double(double)> ConstraintRelaxationStrategy::compute_predicted_objective_reduction_model(const Iterate& current_iterate,
-         const Vector<double>& primal_direction, double step_length) const {
-      // predicted objective reduction: "-∇f(x)^T (αd)"
-      const double directional_derivative = dot(primal_direction, current_iterate.evaluations.objective_gradient);
-      return [=](double objective_multiplier) {
-         return step_length * (-objective_multiplier*directional_derivative) ;
-      };
-   }
-
    void ConstraintRelaxationStrategy::compute_progress_measures(Iterate& current_iterate, Iterate& trial_iterate) {
       if (this->inequality_handling_method->subproblem_definition_changed) {
          DEBUG << "The subproblem definition changed, the globalization strategy is reset and the auxiliary measure is recomputed\n";

uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.hpp

@@ -80,8 +80,6 @@ namespace uno {
       void set_infeasibility_measure(Iterate& iterate) const;
       [[nodiscard]] double compute_predicted_infeasibility_reduction_model(const Iterate& current_iterate, const Vector<double>& primal_direction,
             double step_length) const;
-      [[nodiscard]] std::function<double(double)> compute_predicted_objective_reduction_model(const Iterate& current_iterate,
-            const Vector<double>& primal_direction, double step_length, const SymmetricMatrix<size_t, double>& hessian) const;
       [[nodiscard]] std::function<double(double)> compute_predicted_objective_reduction_model(const Iterate& current_iterate,
             const Vector<double>& primal_direction, double step_length) const;
       void compute_progress_measures(Iterate& current_iterate, Iterate& trial_iterate);

uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp

@@ -207,8 +207,7 @@ namespace uno {
    ProgressMeasures FeasibilityRestoration::compute_predicted_reduction_models(Iterate& current_iterate, const Direction& direction, double step_length) {
       return {
          this->compute_predicted_infeasibility_reduction_model(current_iterate, direction.primals, step_length),
-         this->first_order_predicted_reduction ? this->compute_predicted_objective_reduction_model(current_iterate, direction.primals, step_length) :
-            this->compute_predicted_objective_reduction_model(current_iterate, direction.primals, step_length, this->inequality_handling_method->get_lagrangian_hessian()),
+         this->compute_predicted_objective_reduction_model(current_iterate, direction.primals, step_length),
          this->inequality_handling_method->compute_predicted_auxiliary_reduction_model(this->model, current_iterate, direction.primals, step_length)
       };
    }

uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp

@@ -151,6 +151,10 @@ namespace uno {
          double current_penalty_parameter, WarmstartInformation& warmstart_information) {
       this->l1_relaxed_problem.set_objective_multiplier(current_penalty_parameter);
       this->solve_subproblem(statistics, this->l1_relaxed_problem, current_iterate, current_iterate.multipliers, direction, warmstart_information);
+      if (direction.status == SubproblemStatus::UNBOUNDED_PROBLEM) {
+         throw std::runtime_error("l1Relaxation::solve_l1_relaxed_problem: the subproblem is unbounded, this should not happen. If the subproblem "
+            "has curvature, use regularization. If not, use a trust-region method.\n");
+      }
    }
 
    void l1Relaxation::decrease_parameter_aggressively(Iterate& current_iterate, const Direction& direction) {
@@ -277,8 +281,7 @@ namespace uno {
    ProgressMeasures l1Relaxation::compute_predicted_reduction_models(Iterate& current_iterate, const Direction& direction, double step_length) {
       return {
          this->compute_predicted_infeasibility_reduction_model(current_iterate, direction.primals, step_length),
-         this->first_order_predicted_reduction ? this->compute_predicted_objective_reduction_model(current_iterate, direction.primals, step_length) :
-            this->compute_predicted_objective_reduction_model(current_iterate, direction.primals, step_length, this->inequality_handling_method->get_lagrangian_hessian()),
+         this->compute_predicted_objective_reduction_model(current_iterate, direction.primals, step_length),
          this->inequality_handling_method->compute_predicted_auxiliary_reduction_model(this->model, current_iterate, direction.primals, step_length)
       };
    }

uno/ingredients/hessian_models/ConvexifiedHessian.cpp

@@ -1,34 +1,39 @@
 // Copyright (c) 2024 Charlie Vanaret
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
-#include <stdexcept>
 #include "ConvexifiedHessian.hpp"
+#include "ingredients/constraint_relaxation_strategies/OptimizationProblem.hpp"
 #include "ingredients/hessian_models/UnstableRegularization.hpp"
 #include "ingredients/constraint_relaxation_strategies/OptimizationProblem.hpp"
 #include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
 #include "ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolverFactory.hpp"
-#include "tools/Logger.hpp"
+#include "linear_algebra/SymmetricMatrix.hpp"
 #include "options/Options.hpp"
+#include "tools/Logger.hpp"
 #include "tools/Statistics.hpp"
 
 namespace uno {
    ConvexifiedHessian::ConvexifiedHessian(size_t dimension, size_t maximum_number_nonzeros, const Options& options):
-         HessianModel(dimension, maximum_number_nonzeros, options.get_string("sparse_format"), /* use_regularization = */true),
+         HessianModel(),
          // inertia-based convexification needs a linear solver
          linear_solver(SymmetricIndefiniteLinearSolverFactory::create(dimension, maximum_number_nonzeros, options)),
          regularization_initial_value(options.get_double("regularization_initial_value")),
          regularization_increase_factor(options.get_double("regularization_increase_factor")),
          regularization_failure_threshold(options.get_double("regularization_failure_threshold")) {
    }
 
+   void ConvexifiedHessian::initialize_statistics(Statistics& statistics, const Options& options) const {
+      statistics.add_column("regulariz", Statistics::double_width - 4, options.get_int("statistics_regularization_column_order"));
+   }
+
    void ConvexifiedHessian::evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
-         const Vector<double>& constraint_multipliers) {
+         const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) {
       // evaluate Lagrangian Hessian
-      this->hessian.set_dimension(problem.number_variables);
-      problem.evaluate_lagrangian_hessian(primal_variables, constraint_multipliers, this->hessian);
+      hessian.set_dimension(problem.number_variables);
+      problem.evaluate_lagrangian_hessian(primal_variables, constraint_multipliers, hessian);
       this->evaluation_count++;
       // regularize (only on the original variables) to convexify the problem
-      this->regularize(statistics, this->hessian, problem.get_number_original_variables());
+      this->regularize(statistics, hessian, problem.get_number_original_variables());
    }
 
    // Nocedal and Wright, p51

uno/ingredients/hessian_models/ConvexifiedHessian.hpp

@@ -15,8 +15,9 @@ namespace uno {
    public:
       ConvexifiedHessian(size_t dimension, size_t maximum_number_nonzeros, const Options& options);
 
+      void initialize_statistics(Statistics& statistics, const Options& options) const override;
       void evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
-            const Vector<double>& constraint_multipliers) override;
+            const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) override;
 
    protected:
       std::unique_ptr<DirectSymmetricIndefiniteLinearSolver<size_t, double>> linear_solver; /*!< Solver that computes the inertia */

uno/ingredients/hessian_models/ExactHessian.cpp

@@ -3,19 +3,21 @@
 
 #include "ExactHessian.hpp"
 #include "ingredients/constraint_relaxation_strategies/OptimizationProblem.hpp"
+#include "linear_algebra/SymmetricMatrix.hpp"
 #include "options/Options.hpp"
 
 namespace uno {
    // exact Hessian
-   ExactHessian::ExactHessian(size_t dimension, size_t maximum_number_nonzeros, const Options& options) :
-         HessianModel(dimension, maximum_number_nonzeros, options.get_string("sparse_format"), /* use_regularization = */false) {
+   ExactHessian::ExactHessian(): HessianModel() {
    }
 
+   void ExactHessian::initialize_statistics(Statistics& /*statistics*/, const Options& /*options*/) const { }
+
    void ExactHessian::evaluate(Statistics& /*statistics*/, const OptimizationProblem& problem, const Vector<double>& primal_variables,
-         const Vector<double>& constraint_multipliers) {
+         const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) {
       // evaluate Lagrangian Hessian
-      this->hessian.set_dimension(problem.number_variables);
-      problem.evaluate_lagrangian_hessian(primal_variables, constraint_multipliers, this->hessian);
+      hessian.set_dimension(problem.number_variables);
+      problem.evaluate_lagrangian_hessian(primal_variables, constraint_multipliers, hessian);
       this->evaluation_count++;
    }
 } // namespace

uno/ingredients/hessian_models/ExactHessian.hpp

@@ -4,15 +4,13 @@
 #include "HessianModel.hpp"
 
 namespace uno {
-   // forward declaration
-   class Options;
-
    // exact Hessian
    class ExactHessian : public HessianModel {
    public:
-      ExactHessian(size_t dimension, size_t maximum_number_nonzeros, const Options& options);
+      ExactHessian();
 
+      void initialize_statistics(Statistics& statistics, const Options& options) const override;
       void evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
-            const Vector<double>& constraint_multipliers) override;
+            const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) override;
    };
 } // namespace

uno/ingredients/hessian_models/HessianModel.cpp

@@ -4,9 +4,5 @@
 #include "HessianModel.hpp"
 
 namespace uno {
-   HessianModel::HessianModel(size_t dimension, size_t maximum_number_nonzeros, const std::string& sparse_format, bool use_regularization) :
-         hessian(dimension, maximum_number_nonzeros, use_regularization, sparse_format) {
-   }
-
    HessianModel::~HessianModel() { }
 } // namespace

uno/ingredients/hessian_models/HessianModel.hpp

@@ -4,27 +4,28 @@
 #ifndef UNO_HESSIANMODEL_H
 #define UNO_HESSIANMODEL_H
 
-#include <memory>
-#include <vector>
-#include "linear_algebra/SymmetricMatrix.hpp"
+#include <cstddef>
 
 namespace uno {
    // forward declarations
    class OptimizationProblem;
+   class Options;
    class Statistics;
+   template <typename IndexType, typename ElementType>
+   class SymmetricMatrix;
    template <typename ElementType>
    class Vector;
 
    class HessianModel {
    public:
-      HessianModel(size_t dimension, size_t maximum_number_nonzeros, const std::string& sparse_format, bool use_regularization);
+      HessianModel() = default;
       virtual ~HessianModel();
 
-      SymmetricMatrix<size_t, double> hessian;
       size_t evaluation_count{0};
 
+      virtual void initialize_statistics(Statistics& statistics, const Options& options) const = 0;
       virtual void evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
-            const Vector<double>& constraint_multipliers) = 0;
+            const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) = 0;
    };
 } // namespace
 

uno/ingredients/hessian_models/HessianModelFactory.cpp

@@ -1,6 +1,7 @@
 // Copyright (c) 2018-2024 Charlie Vanaret
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
+#include <stdexcept>
 #include "HessianModelFactory.hpp"
 #include "HessianModel.hpp"
 #include "ConvexifiedHessian.hpp"
@@ -16,11 +17,11 @@ namespace uno {
             return std::make_unique<ConvexifiedHessian>(dimension, maximum_number_nonzeros + dimension, options);
          }
          else {
-            return std::make_unique<ExactHessian>(dimension, maximum_number_nonzeros, options);
+            return std::make_unique<ExactHessian>();
          }
       }
       else if (hessian_model == "zero") {
-         return std::make_unique<ZeroHessian>(dimension, options);
+         return std::make_unique<ZeroHessian>();
       }
       throw std::invalid_argument("Hessian model " + hessian_model + " does not exist");
    }

uno/ingredients/hessian_models/ZeroHessian.cpp

@@ -3,14 +3,15 @@
 
 #include "ZeroHessian.hpp"
 #include "ingredients/constraint_relaxation_strategies/OptimizationProblem.hpp"
+#include "linear_algebra/SymmetricMatrix.hpp"
 #include "options/Options.hpp"
 
 namespace uno {
-   ZeroHessian::ZeroHessian(size_t dimension, const Options& options) :
-         HessianModel(dimension, 0, options.get_string("sparse_format"), /* use_regularization = */false) { }
+   void ZeroHessian::initialize_statistics(Statistics& /*statistics*/, const Options& /*options*/) const { }
 
    void ZeroHessian::evaluate(Statistics& /*statistics*/, const OptimizationProblem& problem, const Vector<double>& /*primal_variables*/,
-         const Vector<double>& /*constraint_multipliers*/) {
-      this->hessian.set_dimension(problem.number_variables);
+         const Vector<double>& /*constraint_multipliers*/, SymmetricMatrix<size_t, double>& hessian) {
+      hessian.set_dimension(problem.number_variables);
+      hessian.reset();
    }
 }

uno/ingredients/hessian_models/ZeroHessian.hpp

@@ -4,15 +4,13 @@
 #include "HessianModel.hpp"
 
 namespace uno {
-   // forward declaration
-   class Options;
-
    // zero Hessian
    class ZeroHessian : public HessianModel {
    public:
-      ZeroHessian(size_t dimension, const Options& options);
+      ZeroHessian() = default;
 
+      void initialize_statistics(Statistics& statistics, const Options& options) const override;
       void evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
-            const Vector<double>& constraint_multipliers) override;
+            const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) override;
    };
 } // namespace

uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp

@@ -16,10 +16,6 @@ namespace uno {
       this->trust_region_radius = new_trust_region_radius;
    }
 
-   const SymmetricMatrix<size_t, double>& InequalityHandlingMethod::get_lagrangian_hessian() const {
-      return this->hessian_model->hessian;
-   }
-
    size_t InequalityHandlingMethod::get_hessian_evaluation_count() const {
       return this->hessian_model->evaluation_count;
    }

uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp

@@ -43,7 +43,7 @@ namespace uno {
       virtual void exit_feasibility_problem(const OptimizationProblem& problem, Iterate& trial_iterate) = 0;
 
       // progress measures
-      [[nodiscard]] const SymmetricMatrix<size_t, double>& get_lagrangian_hessian() const;
+      [[nodiscard]] virtual double hessian_quadratic_product(const Vector<double>& primal_direction) const = 0;
       virtual void set_auxiliary_measure(const Model& model, Iterate& iterate) = 0;
       [[nodiscard]] virtual double compute_predicted_auxiliary_reduction_model(const Model& model, const Iterate& current_iterate,
             const Vector<double>& primal_direction, double step_length) const = 0;

uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp

@@ -23,7 +23,8 @@ namespace uno {
          constraint_jacobian(number_constraints, number_variables) {
    }
 
-   void InequalityConstrainedMethod::initialize_statistics(Statistics& /*statistics*/, const Options& /*options*/) {
+   void InequalityConstrainedMethod::initialize_statistics(Statistics& statistics, const Options& options) {
+      this->hessian_model->initialize_statistics(statistics, options);
    }
 
    void InequalityConstrainedMethod::set_initial_point(const Vector<double>& point) {

uno/ingredients/inequality_handling_methods/inequality_constrained_methods/LPSubproblem.cpp

@@ -14,49 +14,24 @@ namespace uno {
          size_t number_jacobian_nonzeros, const Options& options) :
          InequalityConstrainedMethod("zero", number_variables, number_constraints, 0, false, options),
          solver(LPSolverFactory::create(number_variables, number_constraints,
-               number_objective_gradient_nonzeros, number_jacobian_nonzeros, options)),
-         zero_hessian(SymmetricMatrix<size_t, double>::zero(number_variables)) {
+               number_objective_gradient_nonzeros, number_jacobian_nonzeros, options)) {
    }
 
    LPSubproblem::~LPSubproblem() { }
 
    void LPSubproblem::generate_initial_iterate(const OptimizationProblem& /*problem*/, Iterate& /*initial_iterate*/) {
    }
 
-   void LPSubproblem::evaluate_functions(const OptimizationProblem& problem, Iterate& current_iterate, const WarmstartInformation& warmstart_information) {
-      // objective gradient
-      if (warmstart_information.objective_changed) {
-         problem.evaluate_objective_gradient(current_iterate, this->objective_gradient);
-      }
-      // constraints and constraint Jacobian
-      if (warmstart_information.constraints_changed) {
-         problem.evaluate_constraints(current_iterate, this->constraints);
-         problem.evaluate_constraint_jacobian(current_iterate, this->constraint_jacobian);
-      }
-   }
-
-   void LPSubproblem::solve(Statistics& /*statistics*/, const OptimizationProblem& problem, Iterate& current_iterate,  const Multipliers& current_multipliers,
-         Direction& direction, WarmstartInformation& warmstart_information) {
-      // evaluate the functions at the current iterate
-      this->evaluate_functions(problem, current_iterate, warmstart_information);
-
-      // set bounds of the variable displacements
-      if (warmstart_information.variable_bounds_changed) {
-         this->set_direction_bounds(problem, current_iterate);
-      }
-
-      // set bounds of the linearized constraints
-      if (warmstart_information.constraint_bounds_changed) {
-         this->set_linearized_constraint_bounds(problem, this->constraints);
-      }
-
-      // solve the LP
-      this->solver->solve_LP(problem.number_variables, problem.number_constraints, this->direction_lower_bounds, this->direction_upper_bounds,
-            this->linearized_constraints_lower_bounds, this->linearized_constraints_upper_bounds, this->objective_gradient,
-            this->constraint_jacobian, this->initial_point, direction, warmstart_information);
+   void LPSubproblem::solve(Statistics& /*statistics*/, const OptimizationProblem& problem, Iterate& current_iterate,
+         const Multipliers& current_multipliers, Direction& direction, WarmstartInformation& warmstart_information) {
+      this->solver->solve_LP(problem, current_iterate, this->initial_point, direction, this->trust_region_radius, warmstart_information);
       InequalityConstrainedMethod::compute_dual_displacements(current_multipliers, direction.multipliers);
       this->number_subproblems_solved++;
       // reset the initial point
       this->initial_point.fill(0.);
    }
+
+   double LPSubproblem::hessian_quadratic_product(const Vector<double>& /*primal_direction*/) const {
+      return 0.;
+   }
 } // namespace