diff --git a/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp b/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp index 1e96bf78..f531149a 100644 --- a/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp @@ -72,24 +72,15 @@ namespace uno { } std::function ConstraintRelaxationStrategy::compute_predicted_objective_reduction_model(const Iterate& current_iterate, - const Vector& primal_direction, double step_length, const SymmetricMatrix& hessian) const { + const Vector& 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 ConstraintRelaxationStrategy::compute_predicted_objective_reduction_model(const Iterate& current_iterate, - const Vector& 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"; diff --git a/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.hpp b/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.hpp index 29e4b33c..6fb93ef0 100644 --- a/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.hpp +++ b/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& primal_direction, double step_length) const; - [[nodiscard]] std::function compute_predicted_objective_reduction_model(const Iterate& current_iterate, - const Vector& primal_direction, double step_length, const SymmetricMatrix& hessian) const; [[nodiscard]] std::function compute_predicted_objective_reduction_model(const Iterate& current_iterate, const Vector& primal_direction, double step_length) const; void compute_progress_measures(Iterate& current_iterate, Iterate& trial_iterate); diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp index 85730f0e..7afb32bf 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp +++ b/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) }; } diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp index 4722d4b9..cdcf1b30 100644 --- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp +++ b/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) }; } diff --git a/uno/ingredients/hessian_models/ConvexifiedHessian.cpp b/uno/ingredients/hessian_models/ConvexifiedHessian.cpp index 4964be16..42aefaab 100644 --- a/uno/ingredients/hessian_models/ConvexifiedHessian.cpp +++ b/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 #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& primal_variables, - const Vector& constraint_multipliers) { + const Vector& constraint_multipliers, SymmetricMatrix& 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 diff --git a/uno/ingredients/hessian_models/ConvexifiedHessian.hpp b/uno/ingredients/hessian_models/ConvexifiedHessian.hpp index 9cd2bfe4..c28803e8 100644 --- a/uno/ingredients/hessian_models/ConvexifiedHessian.hpp +++ b/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& primal_variables, - const Vector& constraint_multipliers) override; + const Vector& constraint_multipliers, SymmetricMatrix& hessian) override; protected: std::unique_ptr> linear_solver; /*!< Solver that computes the inertia */ diff --git a/uno/ingredients/hessian_models/ExactHessian.cpp b/uno/ingredients/hessian_models/ExactHessian.cpp index 52a8d58c..500fdb46 100644 --- a/uno/ingredients/hessian_models/ExactHessian.cpp +++ b/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& primal_variables, - const Vector& constraint_multipliers) { + const Vector& constraint_multipliers, SymmetricMatrix& 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 \ No newline at end of file diff --git a/uno/ingredients/hessian_models/ExactHessian.hpp b/uno/ingredients/hessian_models/ExactHessian.hpp index f2edfe5c..b50c5007 100644 --- a/uno/ingredients/hessian_models/ExactHessian.hpp +++ b/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& primal_variables, - const Vector& constraint_multipliers) override; + const Vector& constraint_multipliers, SymmetricMatrix& hessian) override; }; } // namespace \ No newline at end of file diff --git a/uno/ingredients/hessian_models/HessianModel.cpp b/uno/ingredients/hessian_models/HessianModel.cpp index 7aaacb9f..07e6514f 100644 --- a/uno/ingredients/hessian_models/HessianModel.cpp +++ b/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 \ No newline at end of file diff --git a/uno/ingredients/hessian_models/HessianModel.hpp b/uno/ingredients/hessian_models/HessianModel.hpp index 20dc884e..bab56792 100644 --- a/uno/ingredients/hessian_models/HessianModel.hpp +++ b/uno/ingredients/hessian_models/HessianModel.hpp @@ -4,29 +4,28 @@ #ifndef UNO_HESSIANMODEL_H #define UNO_HESSIANMODEL_H -#include -#include -#include "linear_algebra/SymmetricMatrix.hpp" +#include namespace uno { // forward declarations class OptimizationProblem; class Options; class Statistics; + template + class SymmetricMatrix; template 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 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& primal_variables, - const Vector& constraint_multipliers) = 0; + const Vector& constraint_multipliers, SymmetricMatrix& hessian) = 0; }; } // namespace diff --git a/uno/ingredients/hessian_models/HessianModelFactory.cpp b/uno/ingredients/hessian_models/HessianModelFactory.cpp index a7ee8cde..8c8c6451 100644 --- a/uno/ingredients/hessian_models/HessianModelFactory.cpp +++ b/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 #include "HessianModelFactory.hpp" #include "HessianModel.hpp" #include "ConvexifiedHessian.hpp" @@ -16,11 +17,11 @@ namespace uno { return std::make_unique(dimension, maximum_number_nonzeros + dimension, options); } else { - return std::make_unique(dimension, maximum_number_nonzeros, options); + return std::make_unique(); } } else if (hessian_model == "zero") { - return std::make_unique(dimension, options); + return std::make_unique(); } throw std::invalid_argument("Hessian model " + hessian_model + " does not exist"); } diff --git a/uno/ingredients/hessian_models/ZeroHessian.cpp b/uno/ingredients/hessian_models/ZeroHessian.cpp index 8b6d1236..ce9cc528 100644 --- a/uno/ingredients/hessian_models/ZeroHessian.cpp +++ b/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& /*primal_variables*/, - const Vector& /*constraint_multipliers*/) { - this->hessian.set_dimension(problem.number_variables); + const Vector& /*constraint_multipliers*/, SymmetricMatrix& hessian) { + hessian.set_dimension(problem.number_variables); + hessian.reset(); } } \ No newline at end of file diff --git a/uno/ingredients/hessian_models/ZeroHessian.hpp b/uno/ingredients/hessian_models/ZeroHessian.hpp index 08731801..7032a187 100644 --- a/uno/ingredients/hessian_models/ZeroHessian.hpp +++ b/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& primal_variables, - const Vector& constraint_multipliers) override; + const Vector& constraint_multipliers, SymmetricMatrix& hessian) override; }; } // namespace \ No newline at end of file diff --git a/uno/ingredients/subproblems/Subproblem.cpp b/uno/ingredients/subproblems/Subproblem.cpp index 9811242b..068f1087 100644 --- a/uno/ingredients/subproblems/Subproblem.cpp +++ b/uno/ingredients/subproblems/Subproblem.cpp @@ -16,10 +16,6 @@ namespace uno { this->trust_region_radius = new_trust_region_radius; } - const SymmetricMatrix& Subproblem::get_lagrangian_hessian() const { - return this->hessian_model->hessian; - } - size_t Subproblem::get_hessian_evaluation_count() const { return this->hessian_model->evaluation_count; } diff --git a/uno/ingredients/subproblems/Subproblem.hpp b/uno/ingredients/subproblems/Subproblem.hpp index 6c355eff..b855bf04 100644 --- a/uno/ingredients/subproblems/Subproblem.hpp +++ b/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& get_lagrangian_hessian() const; + [[nodiscard]] virtual double hessian_quadratic_product(const Vector& 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& primal_direction, double step_length) const = 0; diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.cpp b/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.cpp index 712d0bef..67ceae10 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.cpp +++ b/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& /*primal_direction*/) const { + return 0.; + } } // namespace \ No newline at end of file diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp b/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp index 81dad7c5..2b39c102 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp +++ b/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& primal_direction) const override; private: // pointer to allow polymorphism diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp b/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp index 51adb257..6e413879 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp +++ b/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& primal_direction) const { + return this->solver->hessian_quadratic_product(primal_direction); + } } // namespace diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp b/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp index 83db8801..441014ed 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp +++ b/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& primal_direction) const override; protected: const bool enforce_linear_constraints_at_initial_iterate; diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.cpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.cpp index a9b77446..3d39fa34 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.cpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.cpp @@ -24,6 +24,7 @@ namespace uno { objective_gradient(2 * number_variables), // original variables + barrier terms constraints(number_constraints), constraint_jacobian(number_constraints, number_variables), + hessian(number_variables, number_hessian_nonzeros, false, "COO"), augmented_system(options.get_string("sparse_format"), number_variables + number_constraints, number_hessian_nonzeros + number_variables /* diagonal barrier terms for bound constraints */ @@ -138,13 +139,13 @@ namespace uno { } // barrier Lagrangian Hessian if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { - this->hessian_model->evaluate(statistics, barrier_problem, current_iterate.primals, current_multipliers.constraints); + this->hessian_model->evaluate(statistics, barrier_problem, current_iterate.primals, current_multipliers.constraints, this->hessian); } // barrier Lagrangian Hessian if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { // original Lagrangian Hessian - this->hessian_model->evaluate(statistics, barrier_problem, current_iterate.primals, current_multipliers.constraints); + this->hessian_model->evaluate(statistics, barrier_problem, current_iterate.primals, current_multipliers.constraints, this->hessian); // diagonal barrier terms (grouped by variable) for (size_t variable_index: Range(barrier_problem.number_variables)) { @@ -157,7 +158,7 @@ namespace uno { const double distance_to_bound = current_iterate.primals[variable_index] - barrier_problem.variable_upper_bound(variable_index); diagonal_barrier_term += current_multipliers.upper_bounds[variable_index] / distance_to_bound; } - this->hessian_model->hessian.insert(diagonal_barrier_term, variable_index, variable_index); + this->hessian.insert(diagonal_barrier_term, variable_index, variable_index); } } } @@ -197,10 +198,14 @@ namespace uno { direction.subproblem_objective = this->evaluate_subproblem_objective(direction); } + double PrimalDualInteriorPointMethod::hessian_quadratic_product(const Vector& /*primal_direction*/) const { + throw std::runtime_error("PrimalDualInteriorPointMethod::hessian_quadratic_product not implemented"); + } + void PrimalDualInteriorPointMethod::assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem, const Multipliers& current_multipliers, WarmstartInformation& warmstart_information) { // assemble, factorize and regularize the augmented matrix - this->augmented_system.assemble_matrix(this->hessian_model->hessian, this->constraint_jacobian, problem.number_variables, problem.number_constraints); + this->augmented_system.assemble_matrix(this->hessian, this->constraint_jacobian, problem.number_variables, problem.number_constraints); DEBUG << "Testing factorization with regularization factors (0, 0)\n"; this->augmented_system.factorize_matrix(*this->linear_solver, warmstart_information); const double dual_regularization_parameter = std::pow(this->barrier_parameter(), this->parameters.regularization_exponent); @@ -345,7 +350,7 @@ namespace uno { double PrimalDualInteriorPointMethod::evaluate_subproblem_objective(const Direction& direction) const { const double linear_term = dot(direction.primals, this->objective_gradient); - const double quadratic_term = this->hessian_model->hessian.quadratic_product(direction.primals, direction.primals) / 2.; + const double quadratic_term = 0.; // TODO this->sol->hessian.quadratic_product(direction.primals, direction.primals) / 2.; return linear_term + quadratic_term; } diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.hpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.hpp index 5d6760b8..9d789d57 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.hpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.hpp @@ -40,6 +40,7 @@ namespace uno { 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& primal_direction) const override; void set_auxiliary_measure(const Model& model, Iterate& iterate) override; [[nodiscard]] double compute_predicted_auxiliary_reduction_model(const Model& model, const Iterate& current_iterate, @@ -51,6 +52,7 @@ namespace uno { SparseVector objective_gradient; /*!< Sparse Jacobian of the objective */ std::vector constraints; /*!< Constraint values (size \f$m)\f$ */ RectangularMatrix constraint_jacobian; /*!< Sparse Jacobian of the constraints */ + SymmetricMatrix hessian; SymmetricIndefiniteLinearSystem augmented_system; const std::unique_ptr> linear_solver; diff --git a/uno/solvers/BQPD/BQPDSolver.cpp b/uno/solvers/BQPD/BQPDSolver.cpp index 69a79cc4..94269373 100644 --- a/uno/solvers/BQPD/BQPDSolver.cpp +++ b/uno/solvers/BQPD/BQPDSolver.cpp @@ -50,6 +50,8 @@ namespace uno { constraint_jacobian(number_constraints, number_variables), bqpd_jacobian(number_jacobian_nonzeros + number_objective_gradient_nonzeros), // Jacobian + objective gradient bqpd_jacobian_sparsity(number_jacobian_nonzeros + number_objective_gradient_nonzeros + number_constraints + 3), + hessian(number_variables, number_hessian_nonzeros, options.get_string("globalization_mechanism") != "TR" || options.get_bool("convexify_QP"), + "CSC"), kmax(problem_type == BQPDProblemType::QP ? options.get_int("BQPD_kmax") : 0), alp(static_cast(this->mlp)), lp(static_cast(this->mlp)), active_set(number_variables + number_constraints), @@ -83,16 +85,20 @@ namespace uno { double trust_region_radius, const WarmstartInformation& warmstart_information) { if (this->print_subproblem) { DEBUG << "QP:\n"; - DEBUG << "Hessian: " << hessian_model.hessian; } this->set_up_subproblem(problem, current_iterate, initial_point, trust_region_radius, warmstart_information); if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { - hessian_model.evaluate(statistics, problem, current_iterate.primals, current_multipliers); - this->save_hessian_to_local_format(hessian_model.hessian); + hessian_model.evaluate(statistics, problem, current_iterate.primals, current_multipliers, this->hessian); + this->save_hessian_to_local_format(this->hessian); } + DEBUG << "Hessian: " << this->hessian; this->solve_subproblem(problem, initial_point, direction, warmstart_information); } + double BQPDSolver::hessian_quadratic_product(const Vector& primal_direction) const { + return this->hessian.quadratic_product(primal_direction, primal_direction); + } + void BQPDSolver::set_up_subproblem(const OptimizationProblem& problem, Iterate& current_iterate, const Vector& initial_point, double trust_region_radius, const WarmstartInformation& warmstart_information) { // initialize wsc_ common block (Hessian & workspace for BQPD) @@ -173,11 +179,13 @@ namespace uno { const int mode_integer = static_cast(mode); // solve the LP/QP + DEBUG2 << "Running BQPD\n"; BQPD(&n, &m, &this->k, &this->kmax, this->bqpd_jacobian.data(), this->bqpd_jacobian_sparsity.data(), direction.primals.data(), this->lower_bounds.data(), this->upper_bounds.data(), &direction.subproblem_objective, &this->fmin, this->gradient_solution.data(), this->residuals.data(), this->w.data(), this->e.data(), this->active_set.data(), this->alp.data(), this->lp.data(), &this->mlp, &this->peq_solution, this->workspace.data(), this->workspace_sparsity.data(), &mode_integer, &this->ifail, this->info.data(), &this->iprint, &this->nout); + DEBUG2 << "Ran BQPD\n"; const BQPDStatus bqpd_status = BQPDSolver::bqpd_status_from_int(this->ifail); direction.status = BQPDSolver::status_from_bqpd_status(bqpd_status); diff --git a/uno/solvers/BQPD/BQPDSolver.hpp b/uno/solvers/BQPD/BQPDSolver.hpp index 445a8150..53052e0d 100644 --- a/uno/solvers/BQPD/BQPDSolver.hpp +++ b/uno/solvers/BQPD/BQPDSolver.hpp @@ -9,14 +9,11 @@ #include "ingredients/subproblems/SubproblemStatus.hpp" #include "linear_algebra/RectangularMatrix.hpp" #include "linear_algebra/SparseVector.hpp" +#include "linear_algebra/SymmetricMatrix.hpp" #include "linear_algebra/Vector.hpp" #include "solvers/QPSolver.hpp" namespace uno { - // forward declarations - template - class SymmetricMatrix; - // see bqpd.f enum class BQPDStatus { OPTIMAL = 0, @@ -55,6 +52,8 @@ namespace uno { const Vector& initial_point, Direction& direction, HessianModel& hessian_model, double trust_region_radius, const WarmstartInformation& warmstart_information) override; + [[nodiscard]] double hessian_quadratic_product(const Vector& primal_direction) const override; + private: const size_t number_hessian_nonzeros; @@ -64,6 +63,7 @@ namespace uno { RectangularMatrix constraint_jacobian; std::vector bqpd_jacobian{}; std::vector bqpd_jacobian_sparsity{}; + SymmetricMatrix hessian; int kmax{0}, mlp{1000}; size_t mxwk0{2000000}, mxiwk0{500000}; diff --git a/uno/solvers/QPSolver.hpp b/uno/solvers/QPSolver.hpp index 7f38d714..dd0e3419 100644 --- a/uno/solvers/QPSolver.hpp +++ b/uno/solvers/QPSolver.hpp @@ -32,6 +32,8 @@ namespace uno { virtual void solve_QP(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Vector& current_multipliers, const Vector& initial_point, Direction& direction, HessianModel& hessian_model, double trust_region_radius, const WarmstartInformation& warmstart_information) = 0; + + [[nodiscard]] virtual double hessian_quadratic_product(const Vector& primal_direction) const = 0; }; } // namespace