(WIP) Moved the Hessian evaluation from HessianModel to subproblem/su…

…bproblem 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->subproblem->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->subproblem->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->subproblem->get_lagrangian_hessian()),
+         this->compute_predicted_objective_reduction_model(current_iterate, direction.primals, step_length),
          this->subproblem->compute_predicted_auxiliary_reduction_model(this->model, current_iterate, direction.primals, step_length)
       };
    }

uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp

@@ -109,7 +109,7 @@ namespace uno {
             Direction feasibility_direction(direction.number_variables, direction.number_constraints);
             this->solve_subproblem(statistics, this->feasibility_problem, current_iterate, current_iterate.feasibility_multipliers, feasibility_direction,
                   warmstart_information);
-            std::swap(direction.multipliers, direction.feasibility_multipliers);
+            std::swap(feasibility_direction.multipliers, feasibility_direction.feasibility_multipliers);
             const double residual_lowest_violation = this->model.constraint_violation(current_iterate.evaluations.constraints +
                   current_iterate.evaluations.constraint_jacobian * feasibility_direction.primals, Norm::L1);
             DEBUG << "Lowest linearized infeasibility mk(dk): " << residual_lowest_violation << '\n';
@@ -277,8 +277,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->subproblem->get_lagrangian_hessian()),
+         this->compute_predicted_objective_reduction_model(current_iterate, direction.primals, step_length),
          this->subproblem->compute_predicted_auxiliary_reduction_model(this->model, current_iterate, direction.primals, step_length)
       };
    }

uno/ingredients/hessian_models/ConvexifiedHessian.cpp

@@ -1,20 +1,19 @@
 // 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/hessian_models/UnstableRegularization.hpp"
+#include "linear_algebra/SymmetricMatrix.hpp"
+#include "options/Options.hpp"
 #include "reformulation/OptimizationProblem.hpp"
 #include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
 #include "solvers/SymmetricIndefiniteLinearSolverFactory.hpp"
 #include "tools/Logger.hpp"
-#include "options/Options.hpp"
-#include "tools/Infinity.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")),
@@ -27,13 +26,13 @@ namespace uno {
    }
 
    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

@@ -17,7 +17,7 @@ namespace uno {
 
       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

@@ -2,22 +2,22 @@
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
 #include "ExactHessian.hpp"
-#include "reformulation/OptimizationProblem.hpp"
+#include "linear_algebra/SymmetricMatrix.hpp"
 #include "options/Options.hpp"
+#include "reformulation/OptimizationProblem.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,16 +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,29 +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

@@ -2,17 +2,16 @@
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
 #include "ZeroHessian.hpp"
+#include "linear_algebra/SymmetricMatrix.hpp"
 #include "reformulation/OptimizationProblem.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,16 +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/subproblems/Subproblem.cpp

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

uno/ingredients/subproblems/Subproblem.hpp

@@ -46,7 +46,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/subproblems/inequality_constrained_methods/LPSubproblem.cpp

@@ -30,4 +30,8 @@ namespace uno {
       // reset the initial point
       this->initial_point.fill(0.);
    }
+
+   double LPSubproblem::hessian_quadratic_product(const Vector<double>& /*primal_direction*/) const {
+      return 0.;
+   }
 } // namespace

uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp

@@ -20,6 +20,7 @@ namespace uno {
       void generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) override;
       void solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,  const Multipliers& current_multipliers,
             Direction& direction, WarmstartInformation& warmstart_information) override;
+      [[nodiscard]] double hessian_quadratic_product(const Vector<double>& primal_direction) const override;
 
    private:
       // pointer to allow polymorphism

uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp

@@ -22,7 +22,7 @@ namespace uno {
          // maximum number of Hessian nonzeros = number nonzeros + possible diagonal inertia correction
          solver(QPSolverFactory::create(number_variables, number_constraints, number_objective_gradient_nonzeros, number_jacobian_nonzeros,
                // if the QP solver is used during preprocessing, we need to allocate the Hessian with at least number_variables elements
-               std::max(this->enforce_linear_constraints_at_initial_iterate ? number_variables : 0, hessian_model->hessian.capacity()),
+               std::max(this->enforce_linear_constraints_at_initial_iterate ? number_variables : 0, number_hessian_nonzeros),
                options)) {
    }
 
@@ -43,4 +43,8 @@ namespace uno {
       // reset the initial point
       this->initial_point.fill(0.);
    }
+
+   double QPSubproblem::hessian_quadratic_product(const Vector<double>& primal_direction) const {
+      return this->solver->hessian_quadratic_product(primal_direction);
+   }
 } // namespace

uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp

@@ -20,6 +20,7 @@ namespace uno {
       void generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) override;
       void solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,  const Multipliers& current_multipliers,
             Direction& direction, WarmstartInformation& warmstart_information) override;
+      [[nodiscard]] double hessian_quadratic_product(const Vector<double>& primal_direction) const override;
 
    protected:
       const bool enforce_linear_constraints_at_initial_iterate;