From 50f28e7c6e12a5fa4accc44d248079df6d11dcab Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 13 Nov 2025 22:16:52 +0100 Subject: [PATCH 1/6] Added class Scaling + EvaluationSpace::compute_constraint_jacobian_norms() + creation of scaling object in constraint relaxation strategies --- uno/Uno.cpp | 1 - .../FeasibilityRestoration.cpp | 6 +++ .../FeasibilityRestoration.hpp | 3 ++ .../NoRelaxation.cpp | 7 ++++ .../NoRelaxation.hpp | 3 ++ .../BQPD/BQPDEvaluationSpace.cpp | 10 +++++ .../BQPD/BQPDEvaluationSpace.hpp | 1 + .../subproblem_solvers/BoxLPSolver.hpp | 3 ++ .../subproblem_solvers/COOEvaluationSpace.cpp | 11 ++++++ .../subproblem_solvers/COOEvaluationSpace.hpp | 2 + .../HiGHS/HiGHSEvaluationSpace.cpp | 11 ++++++ .../HiGHS/HiGHSEvaluationSpace.hpp | 1 + uno/optimization/EvaluationSpace.hpp | 1 + uno/optimization/Scaling.cpp | 38 +++++++++++++++++++ uno/optimization/Scaling.hpp | 30 +++++++++++++++ uno/options/DefaultOptions.cpp | 4 +- uno/options/Options.hpp | 2 +- uno/options/Presets.cpp | 2 +- 18 files changed, 131 insertions(+), 5 deletions(-) create mode 100644 uno/optimization/Scaling.cpp create mode 100644 uno/optimization/Scaling.hpp diff --git a/uno/Uno.cpp b/uno/Uno.cpp index d61b10660..f9363fb4f 100644 --- a/uno/Uno.cpp +++ b/uno/Uno.cpp @@ -2,7 +2,6 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include "Uno.hpp" -#include "ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.hpp" #include "ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategyFactory.hpp" #include "ingredients/globalization_mechanisms/GlobalizationMechanismFactory.hpp" #include "ingredients/globalization_strategies/GlobalizationStrategyFactory.hpp" diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp index cf07320ad..8882b922f 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp @@ -73,6 +73,12 @@ namespace uno { ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->optimality_problem, initial_iterate); this->optimality_globalization_strategy->initialize(statistics, initial_iterate, options); this->feasibility_globalization_strategy.initialize(statistics, initial_iterate, options); + + // optional scaling + if (options.get_bool("use_function_scaling")) { + this->scaling.emplace(initial_iterate, this->optimality_inequality_handling_method->get_evaluation_space(), + options.get_double("function_scaling_threshold")); + } } void FeasibilityRestoration::compute_feasible_direction(Statistics& statistics, Iterate& current_iterate, Direction& direction, diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp index fa420dca3..59171083f 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp +++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp @@ -5,6 +5,7 @@ #define UNO_FEASIBILITYRESTORATION_H #include +#include #include "ConstraintRelaxationStrategy.hpp" #include "relaxed_problems/l1RelaxedProblem.hpp" #include "ingredients/globalization_strategies/GlobalizationStrategy.hpp" @@ -12,6 +13,7 @@ #include "ingredients/globalization_strategies/ProgressMeasures.hpp" #include "ingredients/inertia_correction_strategies/InertiaCorrectionStrategy.hpp" #include "linear_algebra/Vector.hpp" +#include "optimization/Scaling.hpp" namespace uno { // forward declaration @@ -56,6 +58,7 @@ namespace uno { std::unique_ptr feasibility_inequality_handling_method; std::unique_ptr optimality_globalization_strategy; l1MeritFunction feasibility_globalization_strategy; + std::optional scaling{}; // the class maintains multipliers for the other phase (feasibility multipliers if we are in the optimality phase, // and vice versa). These multipliers and those of the iterate are swapped whenever we switch phases. Multipliers other_phase_multipliers; diff --git a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp index 85adcc7ae..83db9f952 100644 --- a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp @@ -9,6 +9,7 @@ #include "optimization/Iterate.hpp" #include "optimization/OptimizationProblem.hpp" #include "optimization/WarmstartInformation.hpp" +#include "options/Options.hpp" #include "symbolic/VectorView.hpp" #include "tools/Logger.hpp" @@ -42,6 +43,12 @@ namespace uno { this->problem.evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient, evaluation_space, initial_iterate); this->compute_primal_dual_residuals(this->problem, initial_iterate); this->globalization_strategy.initialize(statistics, initial_iterate, options); + + // optional scaling + if (options.get_bool("use_function_scaling")) { + this->scaling.emplace(initial_iterate, this->inequality_handling_method->get_evaluation_space(), + options.get_double("function_scaling_threshold")); + } } void NoRelaxation::compute_feasible_direction(Statistics& statistics, Iterate& current_iterate, Direction& direction, diff --git a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.hpp b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.hpp index a37f18b53..6edfc9db5 100644 --- a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.hpp +++ b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.hpp @@ -5,11 +5,13 @@ #define UNO_NORELAXATION_H #include +#include #include "ConstraintRelaxationStrategy.hpp" #include "ingredients/globalization_strategies/l1MeritFunction.hpp" #include "ingredients/hessian_models/HessianModel.hpp" #include "ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp" #include "optimization/OptimizationProblem.hpp" +#include "optimization/Scaling.hpp" namespace uno { class NoRelaxation : public ConstraintRelaxationStrategy { @@ -42,6 +44,7 @@ namespace uno { std::unique_ptr hessian_model; std::unique_ptr> inertia_correction_strategy; l1MeritFunction globalization_strategy; + std::optional scaling{}; }; } // namespace diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp index d78898a81..4a7d65457 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp @@ -76,6 +76,16 @@ namespace uno { } } + void BQPDEvaluationSpace::compute_constraint_jacobian_norms(Vector& row_norms) const { + row_norms.fill(0.); + const size_t number_constraint_jacobian_nonzeros = this->jacobian_row_indices.size(); + for (size_t nonzero_index: Range(number_constraint_jacobian_nonzeros)) { + const size_t constraint_index = static_cast(this->jacobian_row_indices[nonzero_index]); + const double derivative = this->jacobian_values[nonzero_index]; + row_norms[constraint_index] = std::max(row_norms[constraint_index], std::abs(derivative)); + } + } + double BQPDEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const { if (subproblem.has_hessian_operator()) { // linear operator // TODO preallocate diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp index 699aaf39a..f24d338e4 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp @@ -25,6 +25,7 @@ namespace uno { void compute_constraint_jacobian_vector_product(const Vector& vector, Vector& result) const override; void compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const override; + void compute_constraint_jacobian_norms(Vector& row_norms) const override; [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const override; void evaluate_functions(const OptimizationProblem& problem, Iterate& current_iterate, const WarmstartInformation& warmstart_information); diff --git a/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp b/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp index bb4a4c866..f180933fa 100644 --- a/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp +++ b/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp @@ -18,6 +18,9 @@ namespace uno { void compute_constraint_jacobian_vector_product(const Vector& /*vector*/, Vector& /*result*/) const override { } void compute_constraint_jacobian_transposed_vector_product(const Vector& /*vector*/, Vector& /*result*/) const override { } + void compute_constraint_jacobian_norms(Vector& row_norms) const override { + row_norms.fill(0.); + } [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const Vector& /*vector*/) const override { return 0.; } diff --git a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp index e7441d192..1c96fec2c 100644 --- a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp +++ b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp @@ -95,6 +95,17 @@ namespace uno { } } + // compute the inf norm of each row + void COOEvaluationSpace::compute_constraint_jacobian_norms(Vector& row_norms) const { + row_norms.fill(0.); + const size_t offset = this->number_hessian_nonzeros; + for (size_t nonzero_index: Range(this->number_jacobian_nonzeros)) { + const size_t constraint_index = static_cast(this->jacobian_row_indices[nonzero_index]); + const double derivative = this->matrix_values[offset + nonzero_index]; + row_norms[constraint_index] = std::max(row_norms[constraint_index], std::abs(derivative)); + } + } + double COOEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const Vector& /*vector*/) const { return 0.; } diff --git a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp index cb471526d..0f23350f0 100644 --- a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp +++ b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp @@ -6,6 +6,7 @@ #include #include +#include "linear_algebra/Norm.hpp" #include "linear_algebra/Vector.hpp" #include "optimization/EvaluationSpace.hpp" #include "../interfaces/C/uno_int.h" @@ -23,6 +24,7 @@ namespace uno { void compute_constraint_jacobian_vector_product(const Vector& vector, Vector& result) const override; void compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const override; + void compute_constraint_jacobian_norms(Vector& row_norms) const override; [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const override; void set_up_linear_system(Statistics& statistics, const Subproblem& subproblem, DirectSymmetricIndefiniteLinearSolver& linear_solver, diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp index d2dbc723d..f90ed770a 100644 --- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp +++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp @@ -75,6 +75,17 @@ namespace uno { } } + void HiGHSEvaluationSpace::compute_constraint_jacobian_norms(Vector& row_norms) const { + row_norms.fill(0.); + const size_t number_constraint_jacobian_nonzeros = this->jacobian_row_indices.size(); + for (size_t nonzero_index: Range(number_constraint_jacobian_nonzeros)) { + const size_t permuted_nonzero_index = this->jacobian_permutation_vector[nonzero_index]; + const size_t constraint_index = static_cast(this->jacobian_row_indices[permuted_nonzero_index]); + const double derivative = this->model.lp_.a_matrix_.value_[nonzero_index]; + row_norms[constraint_index] = std::max(row_norms[constraint_index], std::abs(derivative)); + } + } + double HiGHSEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const Vector& vector) const { double quadratic_product = 0.; const size_t number_hessian_nonzeros = this->hessian_values.size(); diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp index e2ee04fd0..97a9c64c5 100644 --- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp +++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp @@ -27,6 +27,7 @@ namespace uno { void compute_constraint_jacobian_vector_product(const Vector& vector, Vector& result) const override; void compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const override; + void compute_constraint_jacobian_norms(Vector& row_norms) const override; [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const override; void evaluate_functions(Statistics& statistics, const Subproblem& subproblem, const WarmstartInformation& warmstart_information); diff --git a/uno/optimization/EvaluationSpace.hpp b/uno/optimization/EvaluationSpace.hpp index d42a867b7..92b70b212 100644 --- a/uno/optimization/EvaluationSpace.hpp +++ b/uno/optimization/EvaluationSpace.hpp @@ -25,6 +25,7 @@ namespace uno { virtual void compute_constraint_jacobian_vector_product(const Vector& vector, Vector& result) const = 0; virtual void compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const = 0; + virtual void compute_constraint_jacobian_norms(Vector& row_norms) const = 0; [[nodiscard]] virtual double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const = 0; }; } // namespace diff --git a/uno/optimization/Scaling.cpp b/uno/optimization/Scaling.cpp new file mode 100644 index 000000000..ab8f5b058 --- /dev/null +++ b/uno/optimization/Scaling.cpp @@ -0,0 +1,38 @@ +// Copyright (c) 2018-2025 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "Scaling.hpp" +#include "Iterate.hpp" +#include "linear_algebra/Norm.hpp" +#include "linear_algebra/SparseVector.hpp" +#include "linear_algebra/Vector.hpp" +#include "optimization/EvaluationSpace.hpp" +#include "tools/Logger.hpp" + +namespace uno { + Scaling::Scaling(const Iterate& initial_iterate, const EvaluationSpace& evaluation_space, double gradient_threshold): + gradient_threshold(gradient_threshold), objective_scaling(1.), + constraint_scaling(initial_iterate.evaluations.constraints.size(), 1.) { + // objective + this->objective_scaling = std::min(1., this->gradient_threshold / norm_inf(initial_iterate.evaluations.objective_gradient)); + + // constraints + Vector row_norms(initial_iterate.evaluations.constraints.size()); + evaluation_space.compute_constraint_jacobian_norms(row_norms); + for (size_t constraint_index: Range(this->constraint_scaling.size())) { + this->constraint_scaling[constraint_index] = std::min(1., this->gradient_threshold / row_norms[constraint_index]); + } + + DEBUG << "Objective scaling: " << this->objective_scaling << '\n'; + DEBUG << "Constraint scaling: "; print_vector(DEBUG, this->constraint_scaling); + } + + double Scaling::get_objective_scaling() const { + return this->objective_scaling; + } + + double Scaling::get_constraint_scaling(size_t constraint_index) const { + assert(constraint_index < this->constraint_scaling.size() && "The constraint index is not valid."); + return this->constraint_scaling[constraint_index]; + } +} // namespace \ No newline at end of file diff --git a/uno/optimization/Scaling.hpp b/uno/optimization/Scaling.hpp new file mode 100644 index 000000000..a6853bdfc --- /dev/null +++ b/uno/optimization/Scaling.hpp @@ -0,0 +1,30 @@ +// Copyright (c) 2018-2025 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_SCALING_H +#define UNO_SCALING_H + +#include +#include + +namespace uno { + // forward declarations + class EvaluationSpace; + class Iterate; + + class Scaling { + public: + Scaling(const Iterate& initial_iterate, const EvaluationSpace& evaluation_space, double gradient_threshold); + ~Scaling() = default; + + [[nodiscard]] double get_objective_scaling() const; + [[nodiscard]] double get_constraint_scaling(size_t constraint_index) const; + + protected: + const double gradient_threshold; + double objective_scaling; + std::vector constraint_scaling; + }; +} // namespace + +#endif // UNO_SCALING_H \ No newline at end of file diff --git a/uno/options/DefaultOptions.cpp b/uno/options/DefaultOptions.cpp index f7ddf1d5a..6a45b264b 100644 --- a/uno/options/DefaultOptions.cpp +++ b/uno/options/DefaultOptions.cpp @@ -58,8 +58,8 @@ namespace uno { // Hessian model (exact|zero) options.set_string("hessian_model", "exact"); options.set_string("inertia_correction_strategy", "primal"); - // scale the functions (yes|no) - options.set_bool("scale_functions", false); + // use function scaling based on derivatives at the initial point (yes|no) + options.set_bool("use_function_scaling", false); options.set_double("function_scaling_threshold", 100.); // factor scaling options.set_double("function_scaling_factor", 100.); diff --git a/uno/options/Options.hpp b/uno/options/Options.hpp index 66d62e540..e5bd7d4b6 100644 --- a/uno/options/Options.hpp +++ b/uno/options/Options.hpp @@ -82,7 +82,7 @@ namespace uno { {"globalization_strategy", OptionType::STRING}, {"hessian_model", OptionType::STRING}, {"inertia_correction_strategy", OptionType::STRING}, - {"scale_functions", OptionType::BOOL}, + {"use_function_scaling", OptionType::BOOL}, {"function_scaling_threshold", OptionType::DOUBLE}, {"function_scaling_factor", OptionType::DOUBLE}, {"scale_residuals", OptionType::BOOL}, diff --git a/uno/options/Presets.cpp b/uno/options/Presets.cpp index 72316741e..4a2b8def0 100644 --- a/uno/options/Presets.cpp +++ b/uno/options/Presets.cpp @@ -58,7 +58,7 @@ namespace uno { preset_options.set_double("l1_constraint_violation_coefficient", 1000.); preset_options.set_string("progress_norm", "L1"); preset_options.set_string("residual_norm", "INF"); - preset_options.set_bool("scale_functions", true); + preset_options.set_bool("use_function_scaling", true); preset_options.set_double("primal_tolerance", 1e-8); preset_options.set_double("dual_tolerance", 1e-8); preset_options.set_double("loose_primal_tolerance", 1e-6); From c0adfc5e4e2877191f37f81f68cd9916e0d8705a Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 13 Nov 2025 22:56:36 +0100 Subject: [PATCH 2/6] Pass (optional) scaling to Subproblem creation --- .../FeasibilityRestoration.cpp | 9 ++++++--- .../NoRelaxation.cpp | 5 +++-- .../InequalityHandlingMethod.hpp | 11 +++++++---- .../InequalityConstrainedMethod.cpp | 7 +++---- .../InequalityConstrainedMethod.hpp | 9 +++++---- .../InteriorPointMethod.hpp | 15 ++++++++------- uno/ingredients/subproblem/Subproblem.hpp | 2 ++ 7 files changed, 34 insertions(+), 24 deletions(-) diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp index 8882b922f..deef69a8b 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp @@ -163,7 +163,8 @@ namespace uno { const OptimizationProblem& problem, Iterate& current_iterate, Direction& direction, double trust_region_radius, WarmstartInformation& warmstart_information) { direction.set_dimensions(problem.number_variables, problem.number_constraints); - inequality_handling_method.solve(statistics, current_iterate, direction, trust_region_radius, warmstart_information); + inequality_handling_method.solve(statistics, current_iterate, direction, trust_region_radius, this->scaling, + warmstart_information); direction.norm = norm_inf(view(direction.primals, 0, problem.get_number_original_variables())); DEBUG3 << direction << '\n'; } @@ -206,11 +207,13 @@ namespace uno { // determine acceptability, depending on the current phase if (this->current_phase == Phase::OPTIMALITY) { accept_iterate = this->optimality_inequality_handling_method->is_iterate_acceptable(statistics, - *this->optimality_globalization_strategy, current_iterate, trial_iterate, direction, step_length, user_callbacks); + *this->optimality_globalization_strategy, this->scaling, current_iterate, trial_iterate, direction, + step_length, user_callbacks); } else { accept_iterate = this->feasibility_inequality_handling_method->is_iterate_acceptable(statistics, - this->feasibility_globalization_strategy, current_iterate, trial_iterate, direction, step_length, user_callbacks); + this->feasibility_globalization_strategy, this->scaling, current_iterate, trial_iterate, direction, + step_length, user_callbacks); } trial_iterate.status = this->check_termination(model, trial_iterate); diff --git a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp index 83db9f952..35340a81e 100644 --- a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp @@ -56,7 +56,8 @@ namespace uno { direction.reset(); DEBUG << "Solving the subproblem\n"; direction.set_dimensions(this->problem.number_variables, this->problem.number_constraints); - this->inequality_handling_method->solve(statistics, current_iterate, direction, trust_region_radius, warmstart_information); + this->inequality_handling_method->solve(statistics, current_iterate, direction, trust_region_radius, this->scaling, + warmstart_information); direction.norm = norm_inf(view(direction.primals, 0, this->problem.get_number_original_variables())); DEBUG3 << direction << '\n'; warmstart_information.no_changes(); @@ -75,7 +76,7 @@ namespace uno { Iterate& trial_iterate, const Direction& direction, double step_length, WarmstartInformation& warmstart_information, UserCallbacks& user_callbacks) { const bool accept_iterate = this->inequality_handling_method->is_iterate_acceptable(statistics, this->globalization_strategy, - current_iterate, trial_iterate, direction, step_length, user_callbacks); + this->scaling, current_iterate, trial_iterate, direction, step_length, user_callbacks); trial_iterate.status = this->check_termination(model, trial_iterate); warmstart_information.no_changes(); return accept_iterate; diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp index ac66a1909..ea5a1d50d 100644 --- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp @@ -4,8 +4,10 @@ #ifndef UNO_INEQUALITYHANDLINGMETHOD_H #define UNO_INEQUALITYHANDLINGMETHOD_H +#include #include #include "linear_algebra/Norm.hpp" +#include "optimization/Scaling.hpp" namespace uno { // forward declarations @@ -33,11 +35,12 @@ namespace uno { virtual ~InequalityHandlingMethod() = default; virtual void initialize(const OptimizationProblem& problem, Iterate& current_iterate, - HessianModel& hessian_model, InertiaCorrectionStrategy& inertia_correction_strategy, double trust_region_radius) = 0; + HessianModel& hessian_model, InertiaCorrectionStrategy& inertia_correction_strategy, + double trust_region_radius) = 0; virtual void initialize_statistics(Statistics& statistics, const Options& options) = 0; virtual void generate_initial_iterate(Iterate& initial_iterate) = 0; virtual void solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, double trust_region_radius, - WarmstartInformation& warmstart_information) = 0; + const std::optional& scaling, WarmstartInformation& warmstart_information) = 0; virtual void initialize_feasibility_problem(Iterate& current_iterate) = 0; virtual void set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) = 0; @@ -49,8 +52,8 @@ namespace uno { // progress measures [[nodiscard]] virtual bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy, - Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length, - UserCallbacks& user_callbacks) = 0; + const std::optional& scaling, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, + double step_length, UserCallbacks& user_callbacks) = 0; virtual void postprocess_iterate(Iterate& iterate) = 0; diff --git a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp index 2f79a8a89..594120f3c 100644 --- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp @@ -5,7 +5,6 @@ #include "InequalityConstrainedMethod.hpp" #include "optimization/Iterate.hpp" #include "ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp" -#include "ingredients/hessian_models/HessianModel.hpp" #include "ingredients/subproblem/Subproblem.hpp" #include "ingredients/subproblem_solvers/BoxLPSolverFactory.hpp" #include "ingredients/subproblem_solvers/LPSolverFactory.hpp" @@ -55,7 +54,7 @@ namespace uno { } void InequalityConstrainedMethod::solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, - double trust_region_radius, WarmstartInformation& warmstart_information) { + double trust_region_radius, const std::optional& /*scaling*/, WarmstartInformation& warmstart_information) { // solve the subproblem this->solver->solve(statistics, *this->subproblem, trust_region_radius, this->initial_point, direction, warmstart_information); InequalityConstrainedMethod::compute_dual_displacements(current_iterate.multipliers, direction.multipliers); @@ -99,8 +98,8 @@ namespace uno { } bool InequalityConstrainedMethod::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy, - Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length, - UserCallbacks& user_callbacks) { + const std::optional& /*scaling*/, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, + double step_length, UserCallbacks& user_callbacks) { return InequalityHandlingMethod::is_iterate_acceptable(statistics, globalization_strategy, *this->subproblem, this->get_evaluation_space(), current_iterate, trial_iterate, direction, step_length, user_callbacks); } diff --git a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp index f0c4360cc..ea3fe8d23 100644 --- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp @@ -19,11 +19,12 @@ namespace uno { ~InequalityConstrainedMethod() override = default; void initialize(const OptimizationProblem& problem, Iterate& current_iterate, - HessianModel& hessian_model, InertiaCorrectionStrategy& inertia_correction_strategy, double trust_region_radius) override; + HessianModel& hessian_model, InertiaCorrectionStrategy& inertia_correction_strategy, + double trust_region_radius) override; void initialize_statistics(Statistics& statistics, const Options& options) override; void generate_initial_iterate(Iterate& initial_iterate) override; void solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, double trust_region_radius, - WarmstartInformation& warmstart_information) override; + const std::optional& scaling, WarmstartInformation& warmstart_information) override; void initialize_feasibility_problem(Iterate& current_iterate) override; void set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) override; @@ -35,8 +36,8 @@ namespace uno { // acceptance [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy, - Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length, - UserCallbacks& user_callbacks) override; + const std::optional& scaling, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, + double step_length, UserCallbacks& user_callbacks) override; void postprocess_iterate(Iterate& iterate) override; diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp index 45eb8547c..d5e791b8e 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp @@ -29,7 +29,7 @@ namespace uno { void initialize_statistics(Statistics& statistics, const Options& options) override; void generate_initial_iterate(Iterate& initial_iterate) override; void solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, double trust_region_radius, - WarmstartInformation& warmstart_information) override; + const std::optional& scaling, WarmstartInformation& warmstart_information) override; void initialize_feasibility_problem(Iterate& current_iterate) override; void set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& iterate) override; @@ -41,8 +41,8 @@ namespace uno { // acceptance [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy, - Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length, - UserCallbacks& user_callbacks) override; + const std::optional& scaling, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, + double step_length, UserCallbacks& user_callbacks) override; void postprocess_iterate(Iterate& iterate) override; @@ -128,7 +128,7 @@ namespace uno { template void InteriorPointMethod::solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, - double trust_region_radius, WarmstartInformation& warmstart_information) { + double trust_region_radius, const std::optional& /*scaling*/, WarmstartInformation& warmstart_information) { if (is_finite(trust_region_radius)) { throw std::runtime_error("The interior-point subproblem has a trust region. This is not implemented yet"); } @@ -208,7 +208,8 @@ namespace uno { // where jacobian_coefficient = -1 for p, +1 for n // Note: IPOPT uses a '+' sign because they define the Lagrangian as f(x) + \lambda^T c(x) const double mu = this->barrier_parameter(); - const auto elastic_setting_function = [&](Iterate& iterate, size_t /*constraint_index*/, size_t elastic_index, double jacobian_coefficient) { + const auto elastic_setting_function = [&](Iterate& iterate, size_t /*constraint_index*/, size_t elastic_index, + double jacobian_coefficient) { // precomputations const double constraint_j = 0.; // TODO the constraints are not stored here any more... this->constraints[constraint_index]; const double rho = this->l1_constraint_violation_coefficient; @@ -244,8 +245,8 @@ namespace uno { template bool InteriorPointMethod::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy, - Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length, - UserCallbacks& user_callbacks) { + const std::optional& /*scaling*/, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, + double step_length, UserCallbacks& user_callbacks) { return InequalityHandlingMethod::is_iterate_acceptable(statistics, globalization_strategy, *this->subproblem, this->get_evaluation_space(), current_iterate, trial_iterate, direction, step_length, user_callbacks); } diff --git a/uno/ingredients/subproblem/Subproblem.hpp b/uno/ingredients/subproblem/Subproblem.hpp index 7056096ad..881c9ad15 100644 --- a/uno/ingredients/subproblem/Subproblem.hpp +++ b/uno/ingredients/subproblem/Subproblem.hpp @@ -5,11 +5,13 @@ #define UNO_SUBPROBLEM_H #include +#include #include "optimization/Iterate.hpp" #include "linear_algebra/Matrix.hpp" #include "linear_algebra/MatrixOrder.hpp" #include "optimization/Multipliers.hpp" #include "optimization/OptimizationProblem.hpp" +#include "optimization/Scaling.hpp" #include "symbolic/Range.hpp" #include "symbolic/UnaryNegation.hpp" #include "symbolic/VectorView.hpp" From be2c9464d1fb3c70a18292cdf3cfec7b1498d6d5 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 14 Nov 2025 10:07:30 +0100 Subject: [PATCH 3/6] First draft of using the objective scaling (will fail for ipopt preset) --- interfaces/AMPL/AMPLModel.cpp | 5 +-- .../FeasibilityRestoration.cpp | 8 +++-- .../NoRelaxation.cpp | 6 ++-- .../relaxed_problems/l1RelaxedProblem.cpp | 11 +++--- .../relaxed_problems/l1RelaxedProblem.hpp | 6 ++-- .../hessian_models/ExactHessian.cpp | 6 +++- .../hessian_models/ExactHessian.hpp | 2 +- .../hessian_models/HessianModel.hpp | 5 ++- .../hessian_models/IdentityHessian.cpp | 3 +- .../hessian_models/IdentityHessian.hpp | 3 +- .../hessian_models/ZeroHessian.cpp | 3 +- .../hessian_models/ZeroHessian.hpp | 2 +- .../InequalityHandlingMethod.cpp | 12 ++++--- .../InequalityHandlingMethod.hpp | 8 +++-- .../InequalityConstrainedMethod.cpp | 17 +++++---- .../InequalityConstrainedMethod.hpp | 1 + .../InteriorPointMethod.hpp | 15 +++++--- .../PrimalDualInteriorPointProblem.cpp | 14 ++++---- .../PrimalDualInteriorPointProblem.hpp | 6 ++-- uno/ingredients/subproblem/Subproblem.cpp | 8 ++--- uno/ingredients/subproblem/Subproblem.hpp | 4 +-- .../BQPD/BQPDEvaluationSpace.cpp | 10 +++--- .../BQPD/BQPDEvaluationSpace.hpp | 5 ++- .../subproblem_solvers/BQPD/BQPDSolver.cpp | 35 +++++++++++-------- .../subproblem_solvers/BQPD/BQPDSolver.hpp | 8 ++--- .../subproblem_solvers/BoxLPSolver.cpp | 6 ++-- .../subproblem_solvers/BoxLPSolver.hpp | 4 +-- .../subproblem_solvers/COOEvaluationSpace.cpp | 6 ++-- .../subproblem_solvers/COOEvaluationSpace.hpp | 6 ++-- .../HiGHS/HiGHSEvaluationSpace.cpp | 6 ++-- .../HiGHS/HiGHSEvaluationSpace.hpp | 5 ++- .../subproblem_solvers/HiGHS/HiGHSSolver.cpp | 9 ++--- .../subproblem_solvers/HiGHS/HiGHSSolver.hpp | 6 ++-- .../InequalityConstrainedSolver.hpp | 6 +++- .../subproblem_solvers/LPSolver.hpp | 4 +-- .../subproblem_solvers/MA27/MA27Solver.cpp | 6 ++-- .../subproblem_solvers/MA27/MA27Solver.hpp | 4 +-- .../subproblem_solvers/MA57/MA57Solver.cpp | 6 ++-- .../subproblem_solvers/MA57/MA57Solver.hpp | 4 +-- .../subproblem_solvers/MUMPS/MUMPSSolver.cpp | 6 ++-- .../subproblem_solvers/MUMPS/MUMPSSolver.hpp | 4 +-- .../subproblem_solvers/QPSolver.hpp | 4 +-- .../SymmetricIndefiniteLinearSolver.hpp | 6 ++-- uno/model/Model.hpp | 2 ++ uno/optimization/Iterate.hpp | 2 +- uno/optimization/OptimizationProblem.cpp | 25 +++++++------ uno/optimization/OptimizationProblem.hpp | 18 +++++----- 47 files changed, 205 insertions(+), 143 deletions(-) diff --git a/interfaces/AMPL/AMPLModel.cpp b/interfaces/AMPL/AMPLModel.cpp index 7601a5cb6..ad77d9b3c 100644 --- a/interfaces/AMPL/AMPLModel.cpp +++ b/interfaces/AMPL/AMPLModel.cpp @@ -91,12 +91,13 @@ namespace uno { double AMPLModel::evaluate_objective(const Vector& x) const { fint error_flag = 0; - const double result = this->optimization_sense * (*(this->asl)->p.Objval)(this->asl, 0, const_cast(x.data()), &error_flag); + double objective_value = (*(this->asl)->p.Objval)(this->asl, 0, const_cast(x.data()), &error_flag); if (0 < error_flag) { throw FunctionEvaluationError(); } + objective_value *= this->optimization_sense; ++this->number_model_evaluations.objective; - return result; + return objective_value; } /* diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp index deef69a8b..4dec8cc1f 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp @@ -69,7 +69,7 @@ namespace uno { const auto& evaluation_space = this->optimality_inequality_handling_method->get_evaluation_space(); this->optimality_inequality_handling_method->evaluate_constraint_jacobian(initial_iterate); this->optimality_problem.evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient, - evaluation_space, initial_iterate); + evaluation_space, initial_iterate, this->scaling); ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->optimality_problem, initial_iterate); this->optimality_globalization_strategy->initialize(statistics, initial_iterate, options); this->feasibility_globalization_strategy.initialize(statistics, initial_iterate, options); @@ -79,6 +79,8 @@ namespace uno { this->scaling.emplace(initial_iterate, this->optimality_inequality_handling_method->get_evaluation_space(), options.get_double("function_scaling_threshold")); } + this->optimality_inequality_handling_method->evaluate_progress_measures(initial_iterate, this->scaling); + ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->optimality_problem, initial_iterate); } void FeasibilityRestoration::compute_feasible_direction(Statistics& statistics, Iterate& current_iterate, Direction& direction, @@ -236,13 +238,13 @@ namespace uno { if (this->current_phase == Phase::OPTIMALITY) { this->optimality_problem.evaluate_lagrangian_gradient(iterate.residuals.lagrangian_gradient, - this->optimality_inequality_handling_method->get_evaluation_space(), iterate); + this->optimality_inequality_handling_method->get_evaluation_space(), iterate, this->scaling); ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->optimality_problem, iterate); return ConstraintRelaxationStrategy::check_termination(this->optimality_problem, iterate); } else { this->feasibility_problem.evaluate_lagrangian_gradient(iterate.residuals.lagrangian_gradient, - this->feasibility_inequality_handling_method->get_evaluation_space(), iterate); + this->feasibility_inequality_handling_method->get_evaluation_space(), iterate, this->scaling); ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->feasibility_problem, iterate); return ConstraintRelaxationStrategy::check_termination(this->feasibility_problem, iterate); } diff --git a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp index 35340a81e..193b5b29c 100644 --- a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp @@ -40,7 +40,8 @@ namespace uno { initial_iterate.evaluate_constraints(model); this->inequality_handling_method->evaluate_constraint_jacobian(initial_iterate); const auto& evaluation_space = this->inequality_handling_method->get_evaluation_space(); - this->problem.evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient, evaluation_space, initial_iterate); + this->problem.evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient, evaluation_space, + initial_iterate, this->scaling); this->compute_primal_dual_residuals(this->problem, initial_iterate); this->globalization_strategy.initialize(statistics, initial_iterate, options); @@ -49,6 +50,7 @@ namespace uno { this->scaling.emplace(initial_iterate, this->inequality_handling_method->get_evaluation_space(), options.get_double("function_scaling_threshold")); } + this->inequality_handling_method->evaluate_progress_measures(initial_iterate, this->scaling); } void NoRelaxation::compute_feasible_direction(Statistics& statistics, Iterate& current_iterate, Direction& direction, @@ -87,7 +89,7 @@ namespace uno { iterate.evaluate_constraints(model); const auto& evaluation_space = this->inequality_handling_method->get_evaluation_space(); - this->problem.evaluate_lagrangian_gradient(iterate.residuals.lagrangian_gradient, evaluation_space, iterate); + this->problem.evaluate_lagrangian_gradient(iterate.residuals.lagrangian_gradient, evaluation_space, iterate, this->scaling); ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->problem, iterate); return ConstraintRelaxationStrategy::check_termination(this->problem, iterate); } diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp index af60149a6..c75819896 100644 --- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp @@ -53,13 +53,14 @@ namespace uno { } } - void l1RelaxedProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const { + void l1RelaxedProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient, const std::optional& scaling) const { // scale nabla f(x) by rho if (this->objective_multiplier != 0.) { iterate.evaluate_objective_gradient(this->model); // TODO change this + const double objective_scaling = scaling.has_value() ? scaling->get_objective_scaling() : 1.; for (size_t index: Range(this->model.number_variables)) { - objective_gradient[index] = this->objective_multiplier * iterate.evaluations.objective_gradient[index]; + objective_gradient[index] = this->objective_multiplier * objective_scaling * iterate.evaluations.objective_gradient[index]; } } else { @@ -151,7 +152,7 @@ namespace uno { // Lagrangian gradient split in two parts: objective contribution and constraints' contribution void l1RelaxedProblem::evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, - const EvaluationSpace& evaluation_space, Iterate& iterate) const { + const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional& /*scaling*/) const { lagrangian_gradient.objective_contribution.fill(0.); lagrangian_gradient.constraints_contribution.fill(0.); @@ -191,9 +192,9 @@ namespace uno { } void l1RelaxedProblem::evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector& primal_variables, - const Multipliers& multipliers, double* hessian_values) const { + const Multipliers& multipliers, double* hessian_values, const std::optional& scaling) const { hessian_model.evaluate_hessian(statistics, primal_variables, this->get_objective_multiplier(), - multipliers.constraints, hessian_values); + multipliers.constraints, hessian_values, scaling); // proximal contribution size_t nonzero_index = hessian_model.number_nonzeros(); diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp index b760826fb..957aaf52d 100644 --- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp +++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp @@ -24,7 +24,7 @@ namespace uno { void evaluate_constraints(Iterate& iterate, Vector& constraints) const override; // dense objective gradient - void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const override; + void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient, const std::optional& scaling) const override; // sparsity patterns of Jacobian and Hessian [[nodiscard]] size_t number_jacobian_nonzeros() const override; @@ -38,9 +38,9 @@ namespace uno { // numerical evaluations of Jacobian and Hessian void evaluate_constraint_jacobian(Iterate& iterate, double* jacobian_values) const override; void evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, - const EvaluationSpace& evaluation_space, Iterate& iterate) const override; + const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional& scaling) const override; void evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector& primal_variables, - const Multipliers& multipliers, double* hessian_values) const override; + const Multipliers& multipliers, double* hessian_values, const std::optional& scaling) const override; void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, const Multipliers& multipliers, double* result) const override; diff --git a/uno/ingredients/hessian_models/ExactHessian.cpp b/uno/ingredients/hessian_models/ExactHessian.cpp index 68f1646b2..44ce8c8f2 100644 --- a/uno/ingredients/hessian_models/ExactHessian.cpp +++ b/uno/ingredients/hessian_models/ExactHessian.cpp @@ -34,7 +34,11 @@ namespace uno { } void ExactHessian::evaluate_hessian(Statistics& /*statistics*/, const Vector& primal_variables, - double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values) { + double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values, + const std::optional& scaling) { + if (scaling.has_value()) { + objective_multiplier *= scaling->get_objective_scaling(); + } this->model.evaluate_lagrangian_hessian(primal_variables, objective_multiplier, constraint_multipliers, hessian_values); ++this->evaluation_count; } diff --git a/uno/ingredients/hessian_models/ExactHessian.hpp b/uno/ingredients/hessian_models/ExactHessian.hpp index 9022ce0b9..21b1a1337 100644 --- a/uno/ingredients/hessian_models/ExactHessian.hpp +++ b/uno/ingredients/hessian_models/ExactHessian.hpp @@ -22,7 +22,7 @@ namespace uno { [[nodiscard]] bool is_positive_definite() const override; void evaluate_hessian(Statistics& statistics, const Vector& primal_variables, double objective_multiplier, - const Vector& constraint_multipliers, double* hessian_values) override; + const Vector& constraint_multipliers, double* hessian_values, const std::optional& scaling) override; void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, const Vector& constraint_multipliers, double* result) override; diff --git a/uno/ingredients/hessian_models/HessianModel.hpp b/uno/ingredients/hessian_models/HessianModel.hpp index 202a2d6b7..ab95aa360 100644 --- a/uno/ingredients/hessian_models/HessianModel.hpp +++ b/uno/ingredients/hessian_models/HessianModel.hpp @@ -5,9 +5,11 @@ #define UNO_HESSIANMODEL_H #include +#include #include #include #include "../interfaces/C/uno_int.h" +#include "optimization/Scaling.hpp" namespace uno { // forward declarations @@ -32,7 +34,8 @@ namespace uno { [[nodiscard]] virtual bool is_positive_definite() const = 0; virtual void evaluate_hessian(Statistics& statistics, const Vector& primal_variables, - double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values) = 0; + double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values, + const std::optional& scaling) = 0; virtual void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, const Vector& constraint_multipliers, double* result) = 0; }; diff --git a/uno/ingredients/hessian_models/IdentityHessian.cpp b/uno/ingredients/hessian_models/IdentityHessian.cpp index b7d66228d..a59e0c4b6 100644 --- a/uno/ingredients/hessian_models/IdentityHessian.cpp +++ b/uno/ingredients/hessian_models/IdentityHessian.cpp @@ -38,7 +38,8 @@ namespace uno { } void IdentityHessian::evaluate_hessian(Statistics& /*statistics*/, const Vector& /*primal_variables*/, - double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, double* hessian_values) { + double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, double* hessian_values, + const std::optional& /*scaling*/) { DEBUG << "Setting identity Hessian\n"; for (size_t variable_index: Range(this->number_variables)) { hessian_values[variable_index] = 1.; diff --git a/uno/ingredients/hessian_models/IdentityHessian.hpp b/uno/ingredients/hessian_models/IdentityHessian.hpp index 5fde92d9b..e27a22cc4 100644 --- a/uno/ingredients/hessian_models/IdentityHessian.hpp +++ b/uno/ingredients/hessian_models/IdentityHessian.hpp @@ -19,7 +19,8 @@ namespace uno { [[nodiscard]] bool is_positive_definite() const override; void evaluate_hessian(Statistics& statistics, const Vector& primal_variables, - double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values) override; + double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values, + const std::optional& scaling) override; void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, const Vector& constraint_multipliers, double* result) override; diff --git a/uno/ingredients/hessian_models/ZeroHessian.cpp b/uno/ingredients/hessian_models/ZeroHessian.cpp index ecccc6744..b11c9fc93 100644 --- a/uno/ingredients/hessian_models/ZeroHessian.cpp +++ b/uno/ingredients/hessian_models/ZeroHessian.cpp @@ -33,7 +33,8 @@ namespace uno { } void ZeroHessian::evaluate_hessian(Statistics& /*statistics*/, const Vector& /*primal_variables*/, - double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, double* /*hessian_values*/) { + double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, double* /*hessian_values*/, + const std::optional& /*scaling*/) { // do nothing } diff --git a/uno/ingredients/hessian_models/ZeroHessian.hpp b/uno/ingredients/hessian_models/ZeroHessian.hpp index e8c11858d..ce68afdd2 100644 --- a/uno/ingredients/hessian_models/ZeroHessian.hpp +++ b/uno/ingredients/hessian_models/ZeroHessian.hpp @@ -19,7 +19,7 @@ namespace uno { [[nodiscard]] bool is_positive_definite() const override; void evaluate_hessian(Statistics& statistics, const Vector& primal_variables, double objective_multiplier, - const Vector& constraint_multipliers, double* hessian_values) override; + const Vector& constraint_multipliers, double* hessian_values, const std::optional& scaling) override; void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, const Vector& constraint_multipliers, double* result) override; diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp index f2cca245d..0e216caf1 100644 --- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp @@ -17,15 +17,17 @@ namespace uno { progress_norm(norm_from_string(options.get_string("progress_norm"))) { } - void InequalityHandlingMethod::evaluate_progress_measures(const OptimizationProblem& problem, Iterate& iterate) const { + void InequalityHandlingMethod::evaluate_progress_measures(const OptimizationProblem& problem, Iterate& iterate, + const std::optional& scaling) const { problem.set_infeasibility_measure(iterate, this->progress_norm); - problem.set_objective_measure(iterate); + problem.set_objective_measure(iterate, scaling); problem.set_auxiliary_measure(iterate); } bool InequalityHandlingMethod::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy, - const Subproblem& subproblem, const EvaluationSpace& evaluation_space, Iterate& current_iterate, Iterate& trial_iterate, - const Direction& direction, double step_length, UserCallbacks& user_callbacks) { + const Subproblem& subproblem, const std::optional& scaling, const EvaluationSpace& evaluation_space, + Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length, + UserCallbacks& user_callbacks) { this->postprocess_iterate(trial_iterate); const double objective_multiplier = subproblem.problem.get_objective_multiplier(); @@ -37,7 +39,7 @@ namespace uno { subproblem.problem.set_auxiliary_measure(current_iterate); this->subproblem_definition_changed = false; } - this->evaluate_progress_measures(subproblem.problem, trial_iterate); + this->evaluate_progress_measures(subproblem.problem, trial_iterate, scaling); bool accept_iterate = false; if (direction.norm == 0.) { diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp index ea5a1d50d..4b7d739bd 100644 --- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp @@ -39,6 +39,7 @@ namespace uno { double trust_region_radius) = 0; virtual void initialize_statistics(Statistics& statistics, const Options& options) = 0; virtual void generate_initial_iterate(Iterate& initial_iterate) = 0; + virtual void evaluate_progress_measures(Iterate& iterate, const std::optional& scaling) const = 0; virtual void solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, double trust_region_radius, const std::optional& scaling, WarmstartInformation& warmstart_information) = 0; @@ -68,10 +69,11 @@ namespace uno { // when the parameterization of the subproblem (e.g. penalty or barrier parameter) is updated, signal it bool subproblem_definition_changed{false}; - void evaluate_progress_measures(const OptimizationProblem& problem, Iterate& iterate) const; + void evaluate_progress_measures(const OptimizationProblem& problem, Iterate& iterate, const std::optional& scaling) const; [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy, - const Subproblem& subproblem, const EvaluationSpace& evaluation_space, Iterate& current_iterate, Iterate& trial_iterate, - const Direction& direction, double step_length, UserCallbacks& user_callbacks); + const Subproblem& subproblem, const std::optional& scaling, const EvaluationSpace& evaluation_space, + Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length, + UserCallbacks& user_callbacks); }; } // namespace diff --git a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp index 594120f3c..fa0eb2203 100644 --- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp @@ -49,14 +49,19 @@ namespace uno { // do nothing } - void InequalityConstrainedMethod::generate_initial_iterate(Iterate& initial_iterate) { - this->evaluate_progress_measures(*this->problem, initial_iterate); + void InequalityConstrainedMethod::generate_initial_iterate(Iterate& /*iterate*/) { + // do nothing + } + + void InequalityConstrainedMethod::evaluate_progress_measures(Iterate& iterate, const std::optional& scaling) const { + InequalityHandlingMethod::evaluate_progress_measures(*this->problem, iterate, scaling); } void InequalityConstrainedMethod::solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, - double trust_region_radius, const std::optional& /*scaling*/, WarmstartInformation& warmstart_information) { + double trust_region_radius, const std::optional& scaling, WarmstartInformation& warmstart_information) { // solve the subproblem - this->solver->solve(statistics, *this->subproblem, trust_region_radius, this->initial_point, direction, warmstart_information); + this->solver->solve(statistics, *this->subproblem, trust_region_radius, scaling, this->initial_point, direction, + warmstart_information); InequalityConstrainedMethod::compute_dual_displacements(current_iterate.multipliers, direction.multipliers); ++this->number_subproblems_solved; // reset the initial point @@ -98,10 +103,10 @@ namespace uno { } bool InequalityConstrainedMethod::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy, - const std::optional& /*scaling*/, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, + const std::optional& scaling, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length, UserCallbacks& user_callbacks) { return InequalityHandlingMethod::is_iterate_acceptable(statistics, globalization_strategy, *this->subproblem, - this->get_evaluation_space(), current_iterate, trial_iterate, direction, step_length, user_callbacks); + scaling, this->get_evaluation_space(), current_iterate, trial_iterate, direction, step_length, user_callbacks); } void InequalityConstrainedMethod::postprocess_iterate(Iterate& /*iterate*/) { diff --git a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp index ea3fe8d23..668109009 100644 --- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp @@ -23,6 +23,7 @@ namespace uno { double trust_region_radius) override; void initialize_statistics(Statistics& statistics, const Options& options) override; void generate_initial_iterate(Iterate& initial_iterate) override; + void evaluate_progress_measures(Iterate& iterate, const std::optional& scaling) const override; void solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, double trust_region_radius, const std::optional& scaling, WarmstartInformation& warmstart_information) override; diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp index d5e791b8e..0400e0ba3 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp @@ -28,6 +28,7 @@ namespace uno { InertiaCorrectionStrategy& inertia_correction_strategy, double trust_region_radius) override; void initialize_statistics(Statistics& statistics, const Options& options) override; void generate_initial_iterate(Iterate& initial_iterate) override; + void evaluate_progress_measures(Iterate& iterate, const std::optional& scaling) const override; void solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, double trust_region_radius, const std::optional& scaling, WarmstartInformation& warmstart_information) override; @@ -123,12 +124,16 @@ namespace uno { // resize the initial iterate initial_iterate.set_number_variables(this->barrier_problem->number_variables); this->barrier_problem->generate_initial_iterate(initial_iterate); - this->evaluate_progress_measures(*this->barrier_problem, initial_iterate); + } + + template + void InteriorPointMethod::evaluate_progress_measures(Iterate& iterate, const std::optional& scaling) const { + InequalityHandlingMethod::evaluate_progress_measures(*this->barrier_problem, iterate, scaling); } template void InteriorPointMethod::solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, - double trust_region_radius, const std::optional& /*scaling*/, WarmstartInformation& warmstart_information) { + double trust_region_radius, const std::optional& scaling, WarmstartInformation& warmstart_information) { if (is_finite(trust_region_radius)) { throw std::runtime_error("The interior-point subproblem has a trust region. This is not implemented yet"); } @@ -143,7 +148,7 @@ namespace uno { statistics.set("barrier", this->barrier_parameter()); // compute the primal-dual solution - this->linear_solver->solve_indefinite_system(statistics, *this->subproblem, direction, warmstart_information); + this->linear_solver->solve_indefinite_system(statistics, *this->subproblem, scaling, direction, warmstart_information); ++this->number_subproblems_solved; // check whether the augmented matrix was singular, in which case the subproblem is infeasible @@ -245,10 +250,10 @@ namespace uno { template bool InteriorPointMethod::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy, - const std::optional& /*scaling*/, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, + const std::optional& scaling, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length, UserCallbacks& user_callbacks) { return InequalityHandlingMethod::is_iterate_acceptable(statistics, globalization_strategy, *this->subproblem, - this->get_evaluation_space(), current_iterate, trial_iterate, direction, step_length, user_callbacks); + scaling, this->get_evaluation_space(), current_iterate, trial_iterate, direction, step_length, user_callbacks); } template diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp index 5c71ce0c3..05376e456 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp @@ -69,8 +69,9 @@ namespace uno { this->first_reformulation.evaluate_constraints(iterate, constraints); } - void PrimalDualInteriorPointProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const { - this->first_reformulation.evaluate_objective_gradient(iterate, objective_gradient); + void PrimalDualInteriorPointProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient, + const std::optional& scaling) const { + this->first_reformulation.evaluate_objective_gradient(iterate, objective_gradient, scaling); // barrier terms for (size_t variable_index: Range(this->first_reformulation.number_variables)) { @@ -153,8 +154,8 @@ namespace uno { } void PrimalDualInteriorPointProblem::evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, - const EvaluationSpace& evaluation_space, Iterate& iterate) const { - this->first_reformulation.evaluate_lagrangian_gradient(lagrangian_gradient, evaluation_space, iterate); + const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional& scaling) const { + this->first_reformulation.evaluate_lagrangian_gradient(lagrangian_gradient, evaluation_space, iterate, scaling); // barrier terms for (size_t variable_index: Range(this->first_reformulation.number_variables)) { @@ -179,9 +180,10 @@ namespace uno { } void PrimalDualInteriorPointProblem::evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector& primal_variables, - const Multipliers& multipliers, double* hessian_values) const { + const Multipliers& multipliers, double* hessian_values, const std::optional& scaling) const { // original Lagrangian Hessian - this->first_reformulation.evaluate_lagrangian_hessian(statistics, hessian_model, primal_variables, multipliers, hessian_values); + this->first_reformulation.evaluate_lagrangian_hessian(statistics, hessian_model, primal_variables, multipliers, + hessian_values, scaling); // barrier terms size_t nonzero_index = this->first_reformulation.number_hessian_nonzeros(hessian_model); diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp index 1da605d4b..6c877b679 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp @@ -23,7 +23,7 @@ namespace uno { void evaluate_constraints(Iterate& iterate, Vector& constraints) const override; // dense objective gradient - void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const override; + void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient, const std::optional& scaling) const override; // sparsity patterns of Jacobian and Hessian void compute_constraint_jacobian_sparsity(uno_int* row_indices, uno_int* column_indices, uno_int solver_indexing, @@ -37,9 +37,9 @@ namespace uno { [[nodiscard]] size_t number_hessian_nonzeros(const HessianModel& hessian_model) const override; void evaluate_constraint_jacobian(Iterate& iterate, double* jacobian_values) const override; void evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, - const EvaluationSpace& evaluation_space, Iterate& iterate) const override; + const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional& scaling) const override; void evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector& primal_variables, - const Multipliers& multipliers, double* hessian_values) const override; + const Multipliers& multipliers, double* hessian_values, const std::optional& scaling) const override; void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, const Multipliers& multipliers, double* result) const override; diff --git a/uno/ingredients/subproblem/Subproblem.cpp b/uno/ingredients/subproblem/Subproblem.cpp index cf74d286a..a5b9fb8f9 100644 --- a/uno/ingredients/subproblem/Subproblem.cpp +++ b/uno/ingredients/subproblem/Subproblem.cpp @@ -74,10 +74,10 @@ namespace uno { } } - void Subproblem::evaluate_lagrangian_hessian(Statistics& statistics, double* hessian_values) const { + void Subproblem::evaluate_lagrangian_hessian(Statistics& statistics, double* hessian_values, const std::optional& scaling) const { // evaluate the Lagrangian Hessian of the problem at the current primal-dual point this->problem.evaluate_lagrangian_hessian(statistics, this->hessian_model, this->current_iterate.primals, - this->current_iterate.multipliers, hessian_values); + this->current_iterate.multipliers, hessian_values, scaling); } void Subproblem::regularize_lagrangian_hessian(Statistics& statistics, double* hessian_values) const { @@ -105,10 +105,10 @@ namespace uno { } } - void Subproblem::assemble_augmented_matrix(Statistics& statistics, double* augmented_matrix_values) const { + void Subproblem::assemble_augmented_matrix(Statistics& statistics, double* augmented_matrix_values, const std::optional& scaling) const { // evaluate the Lagrangian Hessian of the problem at the current primal-dual point this->problem.evaluate_lagrangian_hessian(statistics, this->hessian_model, this->current_iterate.primals, - this->current_iterate.multipliers, augmented_matrix_values); + this->current_iterate.multipliers, augmented_matrix_values, scaling); // Jacobian of general constraints this->problem.evaluate_constraint_jacobian(this->current_iterate, augmented_matrix_values + this->number_hessian_nonzeros()); diff --git a/uno/ingredients/subproblem/Subproblem.hpp b/uno/ingredients/subproblem/Subproblem.hpp index 881c9ad15..3bc083fc7 100644 --- a/uno/ingredients/subproblem/Subproblem.hpp +++ b/uno/ingredients/subproblem/Subproblem.hpp @@ -44,12 +44,12 @@ namespace uno { const uno_int* jacobian_row_indices, const uno_int* jacobian_column_indices, uno_int solver_indexing) const; // regularized Hessian - void evaluate_lagrangian_hessian(Statistics& statistics, double* hessian_values) const; + void evaluate_lagrangian_hessian(Statistics& statistics, double* hessian_values, const std::optional& scaling) const; void regularize_lagrangian_hessian(Statistics& statistics, double* hessian_values) const; void compute_hessian_vector_product(const double* x, const double* vector, double* result) const; // augmented system - void assemble_augmented_matrix(Statistics& statistics, double* augmented_matrix_values) const; + void assemble_augmented_matrix(Statistics& statistics, double* augmented_matrix_values, const std::optional& scaling) const; void regularize_augmented_matrix(Statistics& statistics, double* augmented_matrix_values, double dual_regularization_parameter, DirectSymmetricIndefiniteLinearSolver& linear_solver) const; template diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp index 4a7d65457..bc68739fb 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp @@ -115,16 +115,16 @@ namespace uno { } } - void BQPDEvaluationSpace::evaluate_functions(const OptimizationProblem& problem, Iterate& current_iterate, - const WarmstartInformation& warmstart_information) { + void BQPDEvaluationSpace::evaluate_functions(const Subproblem& subproblem, const std::optional& scaling, + Iterate& current_iterate, const WarmstartInformation& warmstart_information) { // evaluate the functions based on warmstart information // gradients is a concatenation of the dense objective gradient and the sparse Jacobian if (warmstart_information.objective_changed) { - problem.evaluate_objective_gradient(current_iterate, this->gradients.data()); + subproblem.problem.evaluate_objective_gradient(current_iterate, this->gradients.data(), scaling); } if (warmstart_information.constraints_changed) { - problem.evaluate_constraints(current_iterate, this->constraints); - this->evaluate_constraint_jacobian(problem, current_iterate); + subproblem.problem.evaluate_constraints(current_iterate, this->constraints); + this->evaluate_constraint_jacobian(subproblem.problem, current_iterate); } if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { this->evaluate_hessian = true; diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp index f24d338e4..15fa30f93 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp @@ -5,8 +5,10 @@ #define UNO_BQPDEVALUATIONSPACE_H #include +#include #include "linear_algebra/Vector.hpp" #include "optimization/EvaluationSpace.hpp" +#include "optimization/Scaling.hpp" #include "../interfaces/C/uno_int.h" namespace uno { @@ -28,7 +30,8 @@ namespace uno { void compute_constraint_jacobian_norms(Vector& row_norms) const override; [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const override; - void evaluate_functions(const OptimizationProblem& problem, Iterate& current_iterate, const WarmstartInformation& warmstart_information); + void evaluate_functions(const Subproblem& subproblem, const std::optional& scaling, Iterate& current_iterate, + const WarmstartInformation& warmstart_information); Vector constraints{}; Vector gradients{}; diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp index 2fa610282..a027332cb 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp @@ -98,17 +98,17 @@ namespace uno { this->mxws = static_cast(this->kmax * (this->kmax + 9) / 2) + 2 * subproblem.number_variables + subproblem.number_constraints /* (required by bqpd.f) */ + 5 * subproblem.number_variables + this->nprof /* (required by sparseL.f) */; - // 6 pointers hidden in lws - constexpr size_t hidden_pointers_size = 6*sizeof(intptr_t); + // 7 pointers hidden in lws + constexpr size_t hidden_pointers_size = 7*sizeof(intptr_t); this->mxlws = hidden_pointers_size + static_cast(this->kmax) /* (required by bqpd.f) */ + 9 * subproblem.number_variables + subproblem.number_constraints /* (required by sparseL.f) */; this->ws.resize(this->mxws); this->lws.resize(this->mxlws); } - void BQPDSolver::solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, + void BQPDSolver::solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional& scaling, const Vector& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) { - this->set_up_subproblem(statistics, subproblem, trust_region_radius, warmstart_information); + this->set_up_subproblem(statistics, subproblem, trust_region_radius, scaling, warmstart_information); if (this->print_subproblem) { this->display_subproblem(subproblem, initial_point); } @@ -122,14 +122,14 @@ namespace uno { // protected member functions void BQPDSolver::set_up_subproblem(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius, - const WarmstartInformation& warmstart_information) { + const std::optional& scaling, const WarmstartInformation& warmstart_information) { // initialize wsc_ common block (Hessian & workspace for BQPD) // setting the common block here ensures that several instances of BQPD can run simultaneously WSC.mxws = static_cast(this->mxws); WSC.mxlws = static_cast(this->mxlws); // evaluate the functions and derivatives - this->evaluation_space.evaluate_functions(subproblem.problem, subproblem.current_iterate, warmstart_information); + this->evaluation_space.evaluate_functions(subproblem, scaling, subproblem.current_iterate, warmstart_information); // variable bounds if (warmstart_information.variable_bounds_changed) { @@ -149,7 +149,7 @@ namespace uno { this->upper_bounds[variable_index] = std::min(BIG, this->upper_bounds[variable_index]); } - this->hide_pointers_in_workspace(statistics, subproblem); + this->hide_pointers_in_workspace(statistics, subproblem, scaling); } void BQPDSolver::display_subproblem(const Subproblem& subproblem, const Vector& initial_point) const { @@ -242,7 +242,8 @@ namespace uno { } // hide pointers to arbitrary objects into this->workspace_sparsity (BQPD's lws) - void BQPDSolver::hide_pointers_in_workspace(Statistics& statistics, const Subproblem& subproblem) { + void BQPDSolver::hide_pointers_in_workspace(Statistics& statistics, const Subproblem& subproblem, + const std::optional& scaling) { WSC.kk = 0; // length of ws that is used by gdotx WSC.ll = 0; // length of lws that is used by gdotx @@ -253,11 +254,13 @@ namespace uno { WSC.ll += sizeof(intptr_t); hide_pointer(2, this->lws.data(), subproblem); WSC.ll += sizeof(intptr_t); - hide_pointer(3, this->lws.data(), this->evaluation_space.hessian_row_indices); + hide_pointer(3, this->lws.data(), scaling); + WSC.ll += sizeof(intptr_t); + hide_pointer(4, this->lws.data(), this->evaluation_space.hessian_row_indices); WSC.ll += sizeof(intptr_t); - hide_pointer(4, this->lws.data(), this->evaluation_space.hessian_column_indices); + hide_pointer(5, this->lws.data(), this->evaluation_space.hessian_column_indices); WSC.ll += sizeof(intptr_t); - hide_pointer(5, this->lws.data(), this->evaluation_space.hessian_values); + hide_pointer(6, this->lws.data(), this->evaluation_space.hessian_values); WSC.ll += sizeof(intptr_t); } @@ -343,12 +346,14 @@ void hessian_vector_product(int* dimension, const double vector[], const double bool* evaluate_hessian = uno::retrieve_pointer(0, lws); uno::Statistics* statistics = uno::retrieve_pointer(1, lws); uno::Subproblem* subproblem = uno::retrieve_pointer(2, lws); - uno::Vector* hessian_row_indices = uno::retrieve_pointer>(3, lws); - uno::Vector* hessian_column_indices = uno::retrieve_pointer>(4, lws); - uno::Vector* hessian_values = uno::retrieve_pointer>(5, lws); + const std::optional* scaling = uno::retrieve_pointer>(3, lws); + uno::Vector* hessian_row_indices = uno::retrieve_pointer>(4, lws); + uno::Vector* hessian_column_indices = uno::retrieve_pointer>(5, lws); + uno::Vector* hessian_values = uno::retrieve_pointer>(6, lws); assert(evaluate_hessian != nullptr); assert(statistics != nullptr); assert(subproblem != nullptr); + assert(scaling != nullptr); assert(hessian_row_indices != nullptr); assert(hessian_column_indices != nullptr); assert(hessian_values != nullptr); @@ -360,7 +365,7 @@ void hessian_vector_product(int* dimension, const double vector[], const double if (subproblem->has_hessian_matrix()) { // if the Hessian has not been evaluated at the current point, evaluate it if (*evaluate_hessian) { - subproblem->evaluate_lagrangian_hessian(*statistics, hessian_values->data()); + subproblem->evaluate_lagrangian_hessian(*statistics, hessian_values->data(), *scaling); subproblem->regularize_lagrangian_hessian(*statistics, hessian_values->data()); *evaluate_hessian = false; } diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.hpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.hpp index d52e95d6d..f6ece1d6f 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.hpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.hpp @@ -45,8 +45,8 @@ namespace uno { void initialize_memory(const Subproblem& subproblem) override; - void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector& initial_point, - Direction& direction, const WarmstartInformation& warmstart_information) override; + void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional& scaling, + const Vector& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) override; [[nodiscard]] EvaluationSpace& get_evaluation_space() override; @@ -74,12 +74,12 @@ namespace uno { const bool print_subproblem; void set_up_subproblem(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius, - const WarmstartInformation& warmstart_information); + const std::optional& scaling, const WarmstartInformation& warmstart_information); void display_subproblem(const Subproblem& subproblem, const Vector& initial_point) const; void solve_subproblem(const Subproblem& subproblem, const Vector& initial_point, Direction& direction, const WarmstartInformation& warmstart_information); [[nodiscard]] static BQPDMode determine_mode(const WarmstartInformation& warmstart_information); - void hide_pointers_in_workspace(Statistics& statistics, const Subproblem& subproblem); + void hide_pointers_in_workspace(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling); void set_multipliers(size_t number_variables, Multipliers& direction_multipliers) const; [[nodiscard]] static BQPDStatus bqpd_status_from_int(int ifail); [[nodiscard]] bool check_sufficient_workspace_size(BQPDStatus bqpd_status); diff --git a/uno/ingredients/subproblem_solvers/BoxLPSolver.cpp b/uno/ingredients/subproblem_solvers/BoxLPSolver.cpp index 8259c0566..fe030801f 100644 --- a/uno/ingredients/subproblem_solvers/BoxLPSolver.cpp +++ b/uno/ingredients/subproblem_solvers/BoxLPSolver.cpp @@ -14,12 +14,14 @@ namespace uno { } void BoxLPSolver::solve(Statistics& /*statistics*/, Subproblem& subproblem, double trust_region_radius, - const Vector& /*initial_point*/, Direction& direction, const WarmstartInformation& /*warmstart_information*/) { + const std::optional& scaling, const Vector& /*initial_point*/, Direction& direction, + const WarmstartInformation& /*warmstart_information*/) { if (0 < subproblem.number_constraints) { throw std::runtime_error("BoxLPSolver cannot solve problems with general constraints"); } // compute the objective gradient - subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->evaluation_space.objective_gradient.data()); + subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->evaluation_space.objective_gradient.data(), + scaling); // compute the variables bounds subproblem.set_variables_bounds(this->variable_lower_bounds, this->variable_upper_bounds, trust_region_radius); diff --git a/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp b/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp index f180933fa..4c352de3b 100644 --- a/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp +++ b/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp @@ -35,8 +35,8 @@ namespace uno { void initialize_memory(const Subproblem& subproblem) override; - void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector& initial_point, - Direction& direction, const WarmstartInformation& warmstart_information) override; + void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional& scaling, + const Vector& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) override; [[nodiscard]] EvaluationSpace& get_evaluation_space() override; diff --git a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp index 1c96fec2c..7bdd1a575 100644 --- a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp +++ b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp @@ -110,11 +110,11 @@ namespace uno { return 0.; } - void COOEvaluationSpace::set_up_linear_system(Statistics& statistics, const Subproblem& subproblem, + void COOEvaluationSpace::set_up_linear_system(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, DirectSymmetricIndefiniteLinearSolver& linear_solver, const WarmstartInformation& warmstart_information) { // evaluate the functions at the current iterate if (warmstart_information.objective_changed) { - subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->objective_gradient.data()); + subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->objective_gradient.data(), scaling); } if (warmstart_information.constraints_changed) { subproblem.problem.evaluate_constraints(subproblem.current_iterate, this->constraints); @@ -128,7 +128,7 @@ namespace uno { this->analysis_performed = true; } // assemble the augmented matrix - subproblem.assemble_augmented_matrix(statistics, this->matrix_values.data()); + subproblem.assemble_augmented_matrix(statistics, this->matrix_values.data(), scaling); // regularize the augmented matrix (this calls the analysis and the factorization) subproblem.regularize_augmented_matrix(statistics, this->matrix_values.data(), subproblem.dual_regularization_factor(), linear_solver); diff --git a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp index 0f23350f0..3a57964b9 100644 --- a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp +++ b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp @@ -5,10 +5,12 @@ #define UNO_COOEVALUATIONSPACE_H #include +#include #include #include "linear_algebra/Norm.hpp" #include "linear_algebra/Vector.hpp" #include "optimization/EvaluationSpace.hpp" +#include "optimization/Scaling.hpp" #include "../interfaces/C/uno_int.h" namespace uno { @@ -27,8 +29,8 @@ namespace uno { void compute_constraint_jacobian_norms(Vector& row_norms) const override; [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const override; - void set_up_linear_system(Statistics& statistics, const Subproblem& subproblem, DirectSymmetricIndefiniteLinearSolver& linear_solver, - const WarmstartInformation& warmstart_information); + void set_up_linear_system(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, + DirectSymmetricIndefiniteLinearSolver& linear_solver, const WarmstartInformation& warmstart_information); Vector objective_gradient{}; /*!< Sparse Jacobian of the objective */ Vector constraints{}; /*!< Constraint values (size \f$m)\f$ */ diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp index f90ed770a..2ab04739f 100644 --- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp +++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp @@ -103,10 +103,10 @@ namespace uno { } void HiGHSEvaluationSpace::evaluate_functions(Statistics& statistics, const Subproblem& subproblem, - const WarmstartInformation& warmstart_information) { + const std::optional& scaling, const WarmstartInformation& warmstart_information) { // evaluate the functions based on warmstart information if (warmstart_information.objective_changed) { - subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->model.lp_.col_cost_.data()); + subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->model.lp_.col_cost_.data(), scaling); } if (warmstart_information.constraints_changed) { subproblem.problem.evaluate_constraints(subproblem.current_iterate, this->constraints); @@ -114,7 +114,7 @@ namespace uno { } // evaluate the Hessian and regularize it if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { - subproblem.evaluate_lagrangian_hessian(statistics, this->hessian_values.data()); + subproblem.evaluate_lagrangian_hessian(statistics, this->hessian_values.data(), scaling); // copy the Hessian with permutation into this->model.hessian_.value_ for (size_t nonzero_index: Range(subproblem.number_regularized_hessian_nonzeros())) { const size_t permuted_nonzero_index = this->hessian_permutation_vector[nonzero_index]; diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp index 97a9c64c5..e4cd8ab24 100644 --- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp +++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp @@ -5,9 +5,11 @@ #define UNO_HIGHSEVALUATIONSPACE_H #include +#include #include "optimization/EvaluationSpace.hpp" #include "Highs.h" #include "linear_algebra/Vector.hpp" +#include "optimization/Scaling.hpp" #include "../interfaces/C/uno_int.h" namespace uno { @@ -30,7 +32,8 @@ namespace uno { void compute_constraint_jacobian_norms(Vector& row_norms) const override; [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const override; - void evaluate_functions(Statistics& statistics, const Subproblem& subproblem, const WarmstartInformation& warmstart_information); + void evaluate_functions(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, + const WarmstartInformation& warmstart_information); HighsModel model; Vector constraints{}; diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp index 1cabd11d3..0783137fd 100644 --- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp +++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp @@ -20,8 +20,9 @@ namespace uno { } void HiGHSSolver::solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, - const Vector& /*initial_point*/, Direction& direction, const WarmstartInformation& warmstart_information) { - this->set_up_subproblem(statistics, subproblem, trust_region_radius, warmstart_information); + const std::optional& scaling, const Vector& /*initial_point*/, Direction& direction, + const WarmstartInformation& warmstart_information) { + this->set_up_subproblem(statistics, subproblem, trust_region_radius, scaling, warmstart_information); this->solve_subproblem(subproblem, direction); } @@ -32,9 +33,9 @@ namespace uno { // protected member functions void HiGHSSolver::set_up_subproblem(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius, - const WarmstartInformation& warmstart_information) { + const std::optional& scaling, const WarmstartInformation& warmstart_information) { // evaluate the functions and derivatives - this->evaluation_space.evaluate_functions(statistics, subproblem, warmstart_information); + this->evaluation_space.evaluate_functions(statistics, subproblem, scaling, warmstart_information); // variable bounds if (warmstart_information.variable_bounds_changed) { diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.hpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.hpp index 5c67a5a35..48109213a 100644 --- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.hpp +++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.hpp @@ -18,8 +18,8 @@ namespace uno { void initialize_memory(const Subproblem& subproblem) override; - void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector& initial_point, - Direction& direction, const WarmstartInformation& warmstart_information) override; + void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional& scaling, + const Vector& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) override; [[nodiscard]] EvaluationSpace& get_evaluation_space() override; @@ -30,7 +30,7 @@ namespace uno { const bool print_subproblem; void set_up_subproblem(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius, - const WarmstartInformation& warmstart_information); + const std::optional& scaling, const WarmstartInformation& warmstart_information); void solve_subproblem(const Subproblem& subproblem, Direction& direction); }; } // namespace diff --git a/uno/ingredients/subproblem_solvers/InequalityConstrainedSolver.hpp b/uno/ingredients/subproblem_solvers/InequalityConstrainedSolver.hpp index 5712efb7f..6d6807a11 100644 --- a/uno/ingredients/subproblem_solvers/InequalityConstrainedSolver.hpp +++ b/uno/ingredients/subproblem_solvers/InequalityConstrainedSolver.hpp @@ -4,6 +4,9 @@ #ifndef UNO_INEQUALITYCONSTRAINEDSOLVER_H #define UNO_INEQUALITYCONSTRAINEDSOLVER_H +#include +#include "optimization/Scaling.hpp" + namespace uno { // forward declarations class Direction; @@ -22,7 +25,8 @@ namespace uno { virtual void initialize_memory(const Subproblem& subproblem) = 0; virtual void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, - const Vector& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) = 0; + const std::optional& scaling, const Vector& initial_point, Direction& direction, + const WarmstartInformation& warmstart_information) = 0; [[nodiscard]] virtual EvaluationSpace& get_evaluation_space() = 0; }; diff --git a/uno/ingredients/subproblem_solvers/LPSolver.hpp b/uno/ingredients/subproblem_solvers/LPSolver.hpp index b2664f47e..80ec85216 100644 --- a/uno/ingredients/subproblem_solvers/LPSolver.hpp +++ b/uno/ingredients/subproblem_solvers/LPSolver.hpp @@ -14,8 +14,8 @@ namespace uno { void initialize_memory(const Subproblem& subproblem) override = 0; - void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector& initial_point, - Direction& direction, const WarmstartInformation& warmstart_information) override = 0; + void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional& scaling, + const Vector& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) override = 0; [[nodiscard]] EvaluationSpace& get_evaluation_space() override = 0; }; diff --git a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp index 6cb3f449a..1774e520e 100644 --- a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp +++ b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp @@ -231,10 +231,10 @@ namespace uno { } } - void MA27Solver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction, - const WarmstartInformation& warmstart_information) { + void MA27Solver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, + Direction& direction, const WarmstartInformation& warmstart_information) { // set up the linear system by evaluating the functions at the current iterate - this->evaluation_space.set_up_linear_system(statistics, subproblem, *this, warmstart_information); + this->evaluation_space.set_up_linear_system(statistics, subproblem, scaling, *this, warmstart_information); // solve the linear system this->solve_indefinite_system(this->evaluation_space.matrix_values, this->evaluation_space.rhs, this->evaluation_space.solution); // assemble the full primal-dual direction diff --git a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp index 3a58f1b91..93cfd182a 100644 --- a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp +++ b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp @@ -45,8 +45,8 @@ namespace uno { void do_symbolic_analysis() override; void do_numerical_factorization(const double* matrix_values) override; void solve_indefinite_system(const Vector& matrix_values, const Vector& rhs, Vector& result) override; - void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction, - const WarmstartInformation& warmstart_information) override; + void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, + Direction& direction, const WarmstartInformation& warmstart_information) override; [[nodiscard]] Inertia get_inertia() const override; [[nodiscard]] size_t number_negative_eigenvalues() const override; diff --git a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp index a5a81ab3e..e74e2cf5b 100644 --- a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp +++ b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp @@ -203,10 +203,10 @@ namespace uno { } } - void MA57Solver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction, - const WarmstartInformation& warmstart_information) { + void MA57Solver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, + Direction& direction, const WarmstartInformation& warmstart_information) { // set up the linear system by evaluating the functions at the current iterate - this->evaluation_space.set_up_linear_system(statistics, subproblem, *this, warmstart_information); + this->evaluation_space.set_up_linear_system(statistics, subproblem, scaling, *this, warmstart_information); // solve the linear system this->solve_indefinite_system(this->evaluation_space.matrix_values, this->evaluation_space.rhs, this->evaluation_space.solution); // assemble the full primal-dual direction diff --git a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp index 13020acc5..921326345 100644 --- a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp +++ b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp @@ -52,8 +52,8 @@ namespace uno { void do_symbolic_analysis() override; void do_numerical_factorization(const double* matrix_values) override; void solve_indefinite_system(const Vector& matrix_values, const Vector& rhs, Vector& result) override; - void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction, - const WarmstartInformation& warmstart_information) override; + void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, + Direction& direction, const WarmstartInformation& warmstart_information) override; [[nodiscard]] Inertia get_inertia() const override; [[nodiscard]] size_t number_negative_eigenvalues() const override; diff --git a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp index dc3c994f9..4341e31b2 100644 --- a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp @@ -94,10 +94,10 @@ namespace uno { dmumps_c(&this->workspace); } - void MUMPSSolver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction, - const WarmstartInformation& warmstart_information) { + void MUMPSSolver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, + Direction& direction, const WarmstartInformation& warmstart_information) { // set up the linear system by evaluating the functions at the current iterate - this->evaluation_space.set_up_linear_system(statistics, subproblem, *this, warmstart_information); + this->evaluation_space.set_up_linear_system(statistics, subproblem, scaling, *this, warmstart_information); // solve the linear system this->solve_indefinite_system(this->evaluation_space.matrix_values, this->evaluation_space.rhs, this->evaluation_space.solution); // assemble the full primal-dual direction diff --git a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp index 06bc93783..ed8ff18a0 100644 --- a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp +++ b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp @@ -21,8 +21,8 @@ namespace uno { void do_symbolic_analysis() override; void do_numerical_factorization(const double* matrix_values) override; void solve_indefinite_system(const Vector& matrix_values, const Vector& rhs, Vector& result) override; - void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction, - const WarmstartInformation& warmstart_information) override; + void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, + Direction& direction, const WarmstartInformation& warmstart_information) override; [[nodiscard]] Inertia get_inertia() const override; [[nodiscard]] size_t number_negative_eigenvalues() const override; diff --git a/uno/ingredients/subproblem_solvers/QPSolver.hpp b/uno/ingredients/subproblem_solvers/QPSolver.hpp index 6e5faf874..6ed02e522 100644 --- a/uno/ingredients/subproblem_solvers/QPSolver.hpp +++ b/uno/ingredients/subproblem_solvers/QPSolver.hpp @@ -14,8 +14,8 @@ namespace uno { void initialize_memory(const Subproblem& subproblem) override = 0; - void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector& initial_point, - Direction& direction, const WarmstartInformation& warmstart_information) override = 0; + void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional& scaling, + const Vector& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) override = 0; }; } // namespace diff --git a/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp b/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp index 20390d524..101dd0216 100644 --- a/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp +++ b/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp @@ -5,6 +5,8 @@ #define UNO_SYMMETRICINDEFINITELINEARSOLVER_H #include +#include +#include "optimization/Scaling.hpp" namespace uno { // forward declarations @@ -27,8 +29,8 @@ namespace uno { virtual void solve_indefinite_system(const Vector& matrix_values, const Vector& rhs, Vector& result) = 0; - virtual void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction, - const WarmstartInformation& warmstart_information) = 0; + virtual void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, + Direction& direction, const WarmstartInformation& warmstart_information) = 0; [[nodiscard]] virtual EvaluationSpace& get_evaluation_space() = 0; }; diff --git a/uno/model/Model.hpp b/uno/model/Model.hpp index 83d562175..f1e240cc3 100644 --- a/uno/model/Model.hpp +++ b/uno/model/Model.hpp @@ -4,11 +4,13 @@ #ifndef UNO_MODEL_H #define UNO_MODEL_H +#include #include #include #include "linear_algebra/MatrixOrder.hpp" #include "linear_algebra/Norm.hpp" #include "optimization/ProblemType.hpp" +#include "optimization/Scaling.hpp" #include "symbolic/VectorExpression.hpp" #include "../interfaces/C/uno_int.h" diff --git a/uno/optimization/Iterate.hpp b/uno/optimization/Iterate.hpp index c0906b5be..9996ee630 100644 --- a/uno/optimization/Iterate.hpp +++ b/uno/optimization/Iterate.hpp @@ -7,8 +7,8 @@ #include "Evaluations.hpp" #include "SolutionStatus.hpp" #include "ingredients/globalization_strategies/ProgressMeasures.hpp" -#include "optimization/Multipliers.hpp" #include "optimization/DualResiduals.hpp" +#include "optimization/Multipliers.hpp" namespace uno { // forward declaration diff --git a/uno/optimization/OptimizationProblem.cpp b/uno/optimization/OptimizationProblem.cpp index 9500e4a50..06df7d29b 100644 --- a/uno/optimization/OptimizationProblem.cpp +++ b/uno/optimization/OptimizationProblem.cpp @@ -30,10 +30,11 @@ namespace uno { constraints = iterate.evaluations.constraints; } - void OptimizationProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const { + void OptimizationProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient, const std::optional& scaling) const { iterate.evaluate_objective_gradient(this->model); + const double objective_scaling = scaling.has_value() ? scaling->get_objective_scaling() : 1.; for (size_t index: Range(this->number_variables)) { - objective_gradient[index] = iterate.evaluations.objective_gradient[index]; + objective_gradient[index] = objective_scaling * iterate.evaluations.objective_gradient[index]; } } @@ -50,12 +51,12 @@ namespace uno { } void OptimizationProblem::compute_constraint_jacobian_sparsity(uno_int *row_indices, uno_int *column_indices, uno_int solver_indexing, - MatrixOrder matrix_order) const { + MatrixOrder matrix_order) const { this->model.compute_constraint_jacobian_sparsity(row_indices, column_indices, solver_indexing, matrix_order); } void OptimizationProblem::compute_hessian_sparsity(const HessianModel& hessian_model, uno_int *row_indices, - uno_int *column_indices, uno_int solver_indexing) const { + uno_int *column_indices, uno_int solver_indexing) const { hessian_model.compute_sparsity(row_indices, column_indices, solver_indexing); } @@ -66,12 +67,12 @@ namespace uno { // Lagrangian gradient ∇f(x_k) - ∇c(x_k) y_k - z_k // split in two parts: objective contribution and constraints' contribution void OptimizationProblem::evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, - const EvaluationSpace& evaluation_space, Iterate& iterate) const { + const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional& scaling) const { lagrangian_gradient.objective_contribution.fill(0.); lagrangian_gradient.constraints_contribution.fill(0.); // ∇f(x_k) - this->evaluate_objective_gradient(iterate, lagrangian_gradient.objective_contribution.data()); + this->evaluate_objective_gradient(iterate, lagrangian_gradient.objective_contribution.data(), scaling); // ∇c(x_k) λ_k evaluation_space.compute_constraint_jacobian_transposed_vector_product(iterate.multipliers.constraints, @@ -86,9 +87,10 @@ namespace uno { } void OptimizationProblem::evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, - const Vector& primal_variables, const Multipliers& multipliers, double* hessian_values) const { + const Vector& primal_variables, const Multipliers& multipliers, double* hessian_values, + const std::optional& scaling) const { hessian_model.evaluate_hessian(statistics, primal_variables, this->get_objective_multiplier(), - multipliers.constraints, hessian_values); + multipliers.constraints, hessian_values, scaling); } void OptimizationProblem::compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, @@ -213,11 +215,12 @@ namespace uno { } // objective measure: scaled objective - void OptimizationProblem::set_objective_measure(Iterate& iterate) const { + void OptimizationProblem::set_objective_measure(Iterate& iterate, const std::optional& scaling) const { iterate.evaluate_objective(this->model); - const double objective = iterate.evaluations.objective; + const double objective_value = iterate.evaluations.objective; + const double objective_scaling = scaling.has_value() ? scaling->get_objective_scaling() : 1.; iterate.progress.objective = [=](double objective_multiplier) { - return objective_multiplier * objective; + return objective_multiplier * objective_scaling * objective_value; }; } diff --git a/uno/optimization/OptimizationProblem.hpp b/uno/optimization/OptimizationProblem.hpp index 3ffaf62ec..4e1c57ce9 100644 --- a/uno/optimization/OptimizationProblem.hpp +++ b/uno/optimization/OptimizationProblem.hpp @@ -37,7 +37,7 @@ namespace uno { virtual void evaluate_constraints(Iterate& iterate, Vector& constraints) const; // dense objective gradient - virtual void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const; + virtual void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient, const std::optional& scaling) const; // sparsity patterns of Jacobian and Hessian [[nodiscard]] virtual size_t number_jacobian_nonzeros() const; @@ -51,9 +51,10 @@ namespace uno { // numerical evaluations of Jacobian and Hessian virtual void evaluate_constraint_jacobian(Iterate& iterate, double* jacobian_values) const; virtual void evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, - const EvaluationSpace& evaluation_space, Iterate& iterate) const; - virtual void evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector& primal_variables, - const Multipliers& multipliers, double* hessian_values) const; + const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional& scaling) const; + virtual void evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, + const Vector& primal_variables, const Multipliers& multipliers, double* hessian_values, + const std::optional& scaling) const; virtual void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, const Multipliers& multipliers, double* result) const; @@ -69,11 +70,12 @@ namespace uno { [[nodiscard]] virtual const Collection& get_inequality_constraints() const; [[nodiscard]] virtual const Collection& get_dual_regularization_constraints() const; - virtual void assemble_primal_dual_direction(const Iterate& current_iterate, const Vector& solution, Direction& direction) const; + virtual void assemble_primal_dual_direction(const Iterate& current_iterate, const Vector& solution, + Direction& direction) const; [[nodiscard]] virtual double dual_regularization_factor() const; - [[nodiscard]] static double stationarity_error(const LagrangianGradient& lagrangian_gradient, double objective_multiplier, - Norm residual_norm); + [[nodiscard]] static double stationarity_error(const LagrangianGradient& lagrangian_gradient, + double objective_multiplier, Norm residual_norm); [[nodiscard]] virtual double complementarity_error(const Vector& primals, const Vector& constraints, const Multipliers& multipliers, double shift_value, Norm residual_norm) const; @@ -82,7 +84,7 @@ namespace uno { // progress measures virtual void set_infeasibility_measure(Iterate& iterate, Norm norm) const; - virtual void set_objective_measure(Iterate& iterate) const; + virtual void set_objective_measure(Iterate& iterate, const std::optional& scaling) const; virtual void set_auxiliary_measure(Iterate& iterate) const; [[nodiscard]] virtual double compute_predicted_auxiliary_reduction(const Iterate& current_iterate, const Vector& primal_direction, double step_length) const; From 5535d741717508065ae893d29963413793ba7eb8 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 14 Nov 2025 10:54:03 +0100 Subject: [PATCH 4/6] Fixed initialization of globalization strategies --- .../FeasibilityRestoration.cpp | 5 ++--- .../constraint_relaxation_strategies/NoRelaxation.cpp | 3 ++- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp index 4dec8cc1f..3b2ed8d45 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp @@ -70,9 +70,6 @@ namespace uno { this->optimality_inequality_handling_method->evaluate_constraint_jacobian(initial_iterate); this->optimality_problem.evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient, evaluation_space, initial_iterate, this->scaling); - ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->optimality_problem, initial_iterate); - this->optimality_globalization_strategy->initialize(statistics, initial_iterate, options); - this->feasibility_globalization_strategy.initialize(statistics, initial_iterate, options); // optional scaling if (options.get_bool("use_function_scaling")) { @@ -81,6 +78,8 @@ namespace uno { } this->optimality_inequality_handling_method->evaluate_progress_measures(initial_iterate, this->scaling); ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->optimality_problem, initial_iterate); + this->optimality_globalization_strategy->initialize(statistics, initial_iterate, options); + this->feasibility_globalization_strategy.initialize(statistics, initial_iterate, options); } void FeasibilityRestoration::compute_feasible_direction(Statistics& statistics, Iterate& current_iterate, Direction& direction, diff --git a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp index 193b5b29c..33afde21a 100644 --- a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp @@ -37,7 +37,6 @@ namespace uno { // initial iterate this->inequality_handling_method->generate_initial_iterate(initial_iterate); initial_iterate.evaluate_objective_gradient(model); - initial_iterate.evaluate_constraints(model); this->inequality_handling_method->evaluate_constraint_jacobian(initial_iterate); const auto& evaluation_space = this->inequality_handling_method->get_evaluation_space(); this->problem.evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient, evaluation_space, @@ -51,6 +50,8 @@ namespace uno { options.get_double("function_scaling_threshold")); } this->inequality_handling_method->evaluate_progress_measures(initial_iterate, this->scaling); + this->compute_primal_dual_residuals(this->problem, initial_iterate); + this->globalization_strategy.initialize(statistics, initial_iterate, options); } void NoRelaxation::compute_feasible_direction(Statistics& statistics, Iterate& current_iterate, Direction& direction, From a46c6d5b96385fb74bc6c06fdcd2de26d3d319dc Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 14 Nov 2025 13:58:23 +0100 Subject: [PATCH 5/6] Fixed Subproblem::compute_predicted_objective_reduction --- .../InequalityHandlingMethod.cpp | 2 +- uno/ingredients/subproblem/Subproblem.cpp | 10 ++++++---- uno/ingredients/subproblem/Subproblem.hpp | 4 ++-- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp index 0e216caf1..dbf30aef0 100644 --- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp @@ -50,7 +50,7 @@ namespace uno { } else { const ProgressMeasures predicted_reductions = subproblem.compute_predicted_reductions(evaluation_space, direction, - step_length, this->progress_norm); + step_length, this->progress_norm, scaling); accept_iterate = globalization_strategy.is_iterate_acceptable(statistics, current_iterate.progress, trial_iterate.progress, predicted_reductions, objective_multiplier); } diff --git a/uno/ingredients/subproblem/Subproblem.cpp b/uno/ingredients/subproblem/Subproblem.cpp index a5b9fb8f9..576fb8eac 100644 --- a/uno/ingredients/subproblem/Subproblem.cpp +++ b/uno/ingredients/subproblem/Subproblem.cpp @@ -241,7 +241,7 @@ namespace uno { } std::function Subproblem::compute_predicted_objective_reduction(const EvaluationSpace& evaluation_space, - const Vector& primal_direction, double step_length) const { + const Vector& primal_direction, double step_length, const std::optional& scaling) const { // predicted objective reduction: "-∇f(x)^T (αd) - α^2/2 d^T H d" const double directional_derivative = dot(primal_direction, this->current_iterate.evaluations.objective_gradient); // if the regularized Hessian is positive definite (as it usually is in line-search methods), we can compute the @@ -249,17 +249,19 @@ namespace uno { const bool is_regularized_hessian_positive_definite = this->hessian_model.is_positive_definite() && this->performs_primal_regularization(); const double quadratic_term = is_regularized_hessian_positive_definite ? 0. : evaluation_space.compute_hessian_quadratic_product(*this, primal_direction); + const double objective_scaling = scaling.has_value() ? scaling->get_objective_scaling() : 1.; return [=](double objective_multiplier) { - return step_length * (-objective_multiplier*directional_derivative) - step_length*step_length/2. * quadratic_term; + return step_length * (-objective_multiplier * objective_scaling * directional_derivative) - + step_length*step_length/2. * quadratic_term; }; } ProgressMeasures Subproblem::compute_predicted_reductions(const EvaluationSpace& evaluation_space, const Direction& direction, - double step_length, Norm norm) const { + double step_length, Norm norm, const std::optional& scaling) const { return { this->compute_predicted_infeasibility_reduction(evaluation_space, this->problem.model, direction.primals, step_length, norm), - this->compute_predicted_objective_reduction(evaluation_space, direction.primals, step_length), + this->compute_predicted_objective_reduction(evaluation_space, direction.primals, step_length, scaling), this->problem.compute_predicted_auxiliary_reduction(this->current_iterate, direction.primals, step_length) }; } diff --git a/uno/ingredients/subproblem/Subproblem.hpp b/uno/ingredients/subproblem/Subproblem.hpp index 3bc083fc7..f7f999bcd 100644 --- a/uno/ingredients/subproblem/Subproblem.hpp +++ b/uno/ingredients/subproblem/Subproblem.hpp @@ -88,9 +88,9 @@ namespace uno { [[nodiscard]] double compute_predicted_infeasibility_reduction(const EvaluationSpace& evaluation_space, const Model& model, const Vector& primal_direction, double step_length, Norm norm) const; [[nodiscard]] std::function compute_predicted_objective_reduction(const EvaluationSpace& evaluation_space, - const Vector& primal_direction, double step_length) const; + const Vector& primal_direction, double step_length, const std::optional& scaling) const; [[nodiscard]] ProgressMeasures compute_predicted_reductions(const EvaluationSpace& evaluation_space, - const Direction& direction, double step_length, Norm norm) const; + const Direction& direction, double step_length, Norm norm, const std::optional& scaling) const; const OptimizationProblem& problem; Iterate& current_iterate; From 0f4a10dedf7422e867ad8a29bb2a86214c10f3fc Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 14 Nov 2025 14:46:18 +0100 Subject: [PATCH 6/6] Added scaling to compute_hessian_vector_product() --- .../relaxed_problems/l1RelaxedProblem.cpp | 5 +++-- .../relaxed_problems/l1RelaxedProblem.hpp | 2 +- uno/ingredients/hessian_models/ExactHessian.cpp | 7 +++++-- uno/ingredients/hessian_models/ExactHessian.hpp | 2 +- uno/ingredients/hessian_models/HessianModel.hpp | 3 ++- uno/ingredients/hessian_models/IdentityHessian.cpp | 3 ++- uno/ingredients/hessian_models/IdentityHessian.hpp | 2 +- uno/ingredients/hessian_models/ZeroHessian.cpp | 3 ++- uno/ingredients/hessian_models/ZeroHessian.hpp | 2 +- .../barrier_problems/PrimalDualInteriorPointProblem.cpp | 4 ++-- .../barrier_problems/PrimalDualInteriorPointProblem.hpp | 2 +- uno/ingredients/subproblem/Subproblem.cpp | 8 +++++--- uno/ingredients/subproblem/Subproblem.hpp | 2 +- .../subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp | 5 +++-- .../subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp | 3 ++- uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp | 2 +- uno/ingredients/subproblem_solvers/BoxLPSolver.hpp | 3 ++- uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp | 3 ++- uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp | 3 ++- .../subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp | 3 ++- .../subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp | 3 ++- uno/optimization/EvaluationSpace.hpp | 6 +++++- uno/optimization/OptimizationProblem.cpp | 5 +++-- uno/optimization/OptimizationProblem.hpp | 2 +- 24 files changed, 52 insertions(+), 31 deletions(-) diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp index c75819896..c8a11b70f 100644 --- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp @@ -209,8 +209,9 @@ namespace uno { } void l1RelaxedProblem::compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, - const Multipliers& multipliers, double* result) const { - hessian_model.compute_hessian_vector_product(x, vector, this->get_objective_multiplier(), multipliers.constraints, result); + const Multipliers& multipliers, double* result, const std::optional& scaling) const { + hessian_model.compute_hessian_vector_product(x, vector, this->get_objective_multiplier(), multipliers.constraints, + result, scaling); // proximal contribution if (this->proximal_center != nullptr && this->proximal_coefficient != 0.) { diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp index 957aaf52d..4e79df31e 100644 --- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp +++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp @@ -42,7 +42,7 @@ namespace uno { void evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector& primal_variables, const Multipliers& multipliers, double* hessian_values, const std::optional& scaling) const override; void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, - const Multipliers& multipliers, double* result) const override; + const Multipliers& multipliers, double* result, const std::optional& scaling) const override; [[nodiscard]] double variable_lower_bound(size_t variable_index) const override; [[nodiscard]] double variable_upper_bound(size_t variable_index) const override; diff --git a/uno/ingredients/hessian_models/ExactHessian.cpp b/uno/ingredients/hessian_models/ExactHessian.cpp index 44ce8c8f2..84f0b15d8 100644 --- a/uno/ingredients/hessian_models/ExactHessian.cpp +++ b/uno/ingredients/hessian_models/ExactHessian.cpp @@ -43,8 +43,11 @@ namespace uno { ++this->evaluation_count; } - void ExactHessian::compute_hessian_vector_product(const double* x, const double* vector, - double objective_multiplier, const Vector& constraint_multipliers, double* result) { + void ExactHessian::compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, + const Vector& constraint_multipliers, double* result, const std::optional& scaling) { + if (scaling.has_value()) { + objective_multiplier *= scaling->get_objective_scaling(); + } this->model.compute_hessian_vector_product(x, vector, objective_multiplier, constraint_multipliers, result); ++this->evaluation_count; } diff --git a/uno/ingredients/hessian_models/ExactHessian.hpp b/uno/ingredients/hessian_models/ExactHessian.hpp index 21b1a1337..9128668c3 100644 --- a/uno/ingredients/hessian_models/ExactHessian.hpp +++ b/uno/ingredients/hessian_models/ExactHessian.hpp @@ -24,7 +24,7 @@ namespace uno { void evaluate_hessian(Statistics& statistics, const Vector& primal_variables, double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values, const std::optional& scaling) override; void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& constraint_multipliers, double* result) override; + const Vector& constraint_multipliers, double* result, const std::optional& scaling) override; protected: const Model& model; diff --git a/uno/ingredients/hessian_models/HessianModel.hpp b/uno/ingredients/hessian_models/HessianModel.hpp index ab95aa360..a1f33e0cb 100644 --- a/uno/ingredients/hessian_models/HessianModel.hpp +++ b/uno/ingredients/hessian_models/HessianModel.hpp @@ -37,7 +37,8 @@ namespace uno { double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values, const std::optional& scaling) = 0; virtual void compute_hessian_vector_product(const double* x, const double* vector, - double objective_multiplier, const Vector& constraint_multipliers, double* result) = 0; + double objective_multiplier, const Vector& constraint_multipliers, double* result, + const std::optional& scaling) = 0; }; } // namespace diff --git a/uno/ingredients/hessian_models/IdentityHessian.cpp b/uno/ingredients/hessian_models/IdentityHessian.cpp index a59e0c4b6..22af946a7 100644 --- a/uno/ingredients/hessian_models/IdentityHessian.cpp +++ b/uno/ingredients/hessian_models/IdentityHessian.cpp @@ -47,7 +47,8 @@ namespace uno { } void IdentityHessian::compute_hessian_vector_product(const double* /*x*/, const double* vector, - double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, double* result) { + double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, double* result, + const std::optional& /*scaling*/) { for (size_t variable_index: Range(this->number_variables)) { result[variable_index] = vector[variable_index]; } diff --git a/uno/ingredients/hessian_models/IdentityHessian.hpp b/uno/ingredients/hessian_models/IdentityHessian.hpp index e27a22cc4..7332e28e5 100644 --- a/uno/ingredients/hessian_models/IdentityHessian.hpp +++ b/uno/ingredients/hessian_models/IdentityHessian.hpp @@ -22,7 +22,7 @@ namespace uno { double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values, const std::optional& scaling) override; void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& constraint_multipliers, double* result) override; + const Vector& constraint_multipliers, double* result, const std::optional& scaling) override; protected: const size_t number_variables; diff --git a/uno/ingredients/hessian_models/ZeroHessian.cpp b/uno/ingredients/hessian_models/ZeroHessian.cpp index b11c9fc93..299aa0471 100644 --- a/uno/ingredients/hessian_models/ZeroHessian.cpp +++ b/uno/ingredients/hessian_models/ZeroHessian.cpp @@ -39,7 +39,8 @@ namespace uno { } void ZeroHessian::compute_hessian_vector_product(const double* /*x*/, const double* /*vector*/, - double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, double* result) { + double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, double* result, + const std::optional& /*scaling*/) { for (size_t variable_index: Range(this->number_variables)) { result[variable_index] = 0.; } diff --git a/uno/ingredients/hessian_models/ZeroHessian.hpp b/uno/ingredients/hessian_models/ZeroHessian.hpp index ce68afdd2..7ce996973 100644 --- a/uno/ingredients/hessian_models/ZeroHessian.hpp +++ b/uno/ingredients/hessian_models/ZeroHessian.hpp @@ -21,7 +21,7 @@ namespace uno { void evaluate_hessian(Statistics& statistics, const Vector& primal_variables, double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values, const std::optional& scaling) override; void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& constraint_multipliers, double* result) override; + const Vector& constraint_multipliers, double* result, const std::optional& scaling) override; protected: const size_t number_variables; diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp index 05376e456..d24f8414f 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp @@ -207,9 +207,9 @@ namespace uno { } void PrimalDualInteriorPointProblem::compute_hessian_vector_product(HessianModel& hessian_model, const double* x, - const double* vector, const Multipliers& multipliers, double* result) const { + const double* vector, const Multipliers& multipliers, double* result, const std::optional& scaling) const { // original Lagrangian Hessian - this->first_reformulation.compute_hessian_vector_product(hessian_model, x, vector, multipliers, result); + this->first_reformulation.compute_hessian_vector_product(hessian_model, x, vector, multipliers, result, scaling); // barrier terms for (size_t variable_index: Range(this->first_reformulation.number_variables)) { diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp index 6c877b679..5f3add3de 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp @@ -41,7 +41,7 @@ namespace uno { void evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector& primal_variables, const Multipliers& multipliers, double* hessian_values, const std::optional& scaling) const override; void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, - const Multipliers& multipliers, double* result) const override; + const Multipliers& multipliers, double* result, const std::optional& scaling) const override; [[nodiscard]] double variable_lower_bound(size_t variable_index) const override; [[nodiscard]] double variable_upper_bound(size_t variable_index) const override; diff --git a/uno/ingredients/subproblem/Subproblem.cpp b/uno/ingredients/subproblem/Subproblem.cpp index 576fb8eac..bdb0b79b4 100644 --- a/uno/ingredients/subproblem/Subproblem.cpp +++ b/uno/ingredients/subproblem/Subproblem.cpp @@ -92,9 +92,11 @@ namespace uno { } } - void Subproblem::compute_hessian_vector_product(const double* x, const double* vector, double* result) const { + void Subproblem::compute_hessian_vector_product(const double* x, const double* vector, double* result, + const std::optional& scaling) const { // unregularized Hessian-vector product - this->problem.compute_hessian_vector_product(this->hessian_model, x, vector, this->current_iterate.multipliers, result); + this->problem.compute_hessian_vector_product(this->hessian_model, x, vector, this->current_iterate.multipliers, + result, scaling); // contribution of the regularization strategy const double regularization_factor = this->inertia_correction_strategy.get_primal_regularization_factor(); @@ -248,7 +250,7 @@ namespace uno { // predicted reduction with only first-order information (the directional derivative) const bool is_regularized_hessian_positive_definite = this->hessian_model.is_positive_definite() && this->performs_primal_regularization(); const double quadratic_term = is_regularized_hessian_positive_definite ? 0. : - evaluation_space.compute_hessian_quadratic_product(*this, primal_direction); + evaluation_space.compute_hessian_quadratic_product(*this, scaling, primal_direction); const double objective_scaling = scaling.has_value() ? scaling->get_objective_scaling() : 1.; return [=](double objective_multiplier) { return step_length * (-objective_multiplier * objective_scaling * directional_derivative) - diff --git a/uno/ingredients/subproblem/Subproblem.hpp b/uno/ingredients/subproblem/Subproblem.hpp index f7f999bcd..3a1d8ca2d 100644 --- a/uno/ingredients/subproblem/Subproblem.hpp +++ b/uno/ingredients/subproblem/Subproblem.hpp @@ -46,7 +46,7 @@ namespace uno { // regularized Hessian void evaluate_lagrangian_hessian(Statistics& statistics, double* hessian_values, const std::optional& scaling) const; void regularize_lagrangian_hessian(Statistics& statistics, double* hessian_values) const; - void compute_hessian_vector_product(const double* x, const double* vector, double* result) const; + void compute_hessian_vector_product(const double* x, const double* vector, double* result, const std::optional& scaling) const; // augmented system void assemble_augmented_matrix(Statistics& statistics, double* augmented_matrix_values, const std::optional& scaling) const; diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp index bc68739fb..ed4d4b691 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp @@ -86,12 +86,13 @@ namespace uno { } } - double BQPDEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const { + double BQPDEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& subproblem, const std::optional& scaling, + const Vector& vector) const { if (subproblem.has_hessian_operator()) { // linear operator // TODO preallocate Vector result(subproblem.number_variables); // compute Hv - subproblem.compute_hessian_vector_product(subproblem.current_iterate.primals.data(), vector.data(), result.data()); + subproblem.compute_hessian_vector_product(subproblem.current_iterate.primals.data(), vector.data(), result.data(), scaling); // compute the dot product return dot(vector, result); } diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp index 15fa30f93..84f91f0c8 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp @@ -28,7 +28,8 @@ namespace uno { void compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const override; void compute_constraint_jacobian_norms(Vector& row_norms) const override; - [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const override; + [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const std::optional& scaling, + const Vector& vector) const override; void evaluate_functions(const Subproblem& subproblem, const std::optional& scaling, Iterate& current_iterate, const WarmstartInformation& warmstart_information); diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp index a027332cb..1b9559944 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp @@ -386,6 +386,6 @@ void hessian_vector_product(int* dimension, const double vector[], const double } // otherwise, try to perform a Hessian-vector product if possible else if (subproblem->has_hessian_operator()) { - subproblem->compute_hessian_vector_product(subproblem->current_iterate.primals.data(), vector, result); + subproblem->compute_hessian_vector_product(subproblem->current_iterate.primals.data(), vector, result, *scaling); } } \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp b/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp index 4c352de3b..4f7fae960 100644 --- a/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp +++ b/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp @@ -21,7 +21,8 @@ namespace uno { void compute_constraint_jacobian_norms(Vector& row_norms) const override { row_norms.fill(0.); } - [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const Vector& /*vector*/) const override { + [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const std::optional& /*scaling*/, + const Vector& /*vector*/) const override { return 0.; } diff --git a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp index 7bdd1a575..12f7e7d9c 100644 --- a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp +++ b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp @@ -106,7 +106,8 @@ namespace uno { } } - double COOEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const Vector& /*vector*/) const { + double COOEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const std::optional& /*scaling*/, + const Vector& /*vector*/) const { return 0.; } diff --git a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp index 3a57964b9..cf031de9a 100644 --- a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp +++ b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp @@ -27,7 +27,8 @@ namespace uno { void compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const override; void compute_constraint_jacobian_norms(Vector& row_norms) const override; - [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const override; + [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const std::optional& scaling, + const Vector& vector) const override; void set_up_linear_system(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, DirectSymmetricIndefiniteLinearSolver& linear_solver, const WarmstartInformation& warmstart_information); diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp index 2ab04739f..cd0e70e4a 100644 --- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp +++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp @@ -86,7 +86,8 @@ namespace uno { } } - double HiGHSEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const Vector& vector) const { + double HiGHSEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const std::optional& /*scaling*/, + const Vector& vector) const { double quadratic_product = 0.; const size_t number_hessian_nonzeros = this->hessian_values.size(); for (size_t nonzero_index: Range(number_hessian_nonzeros)) { diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp index e4cd8ab24..d599e41a7 100644 --- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp +++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp @@ -30,7 +30,8 @@ namespace uno { void compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const override; void compute_constraint_jacobian_norms(Vector& row_norms) const override; - [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const override; + [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const std::optional& scaling, + const Vector& vector) const override; void evaluate_functions(Statistics& statistics, const Subproblem& subproblem, const std::optional& scaling, const WarmstartInformation& warmstart_information); diff --git a/uno/optimization/EvaluationSpace.hpp b/uno/optimization/EvaluationSpace.hpp index 92b70b212..4afa88198 100644 --- a/uno/optimization/EvaluationSpace.hpp +++ b/uno/optimization/EvaluationSpace.hpp @@ -4,6 +4,9 @@ #ifndef UNO_EVALUATIONSPACE_H #define UNO_EVALUATIONSPACE_H +#include +#include "optimization/Scaling.hpp" + namespace uno { // forward declaration template @@ -26,7 +29,8 @@ namespace uno { virtual void compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const = 0; virtual void compute_constraint_jacobian_norms(Vector& row_norms) const = 0; - [[nodiscard]] virtual double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const = 0; + [[nodiscard]] virtual double compute_hessian_quadratic_product(const Subproblem& subproblem, const std::optional& scaling, + const Vector& vector) const = 0; }; } // namespace diff --git a/uno/optimization/OptimizationProblem.cpp b/uno/optimization/OptimizationProblem.cpp index 06df7d29b..4be1c23ce 100644 --- a/uno/optimization/OptimizationProblem.cpp +++ b/uno/optimization/OptimizationProblem.cpp @@ -94,8 +94,9 @@ namespace uno { } void OptimizationProblem::compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, - const Multipliers& multipliers, double* result) const { - hessian_model.compute_hessian_vector_product(x, vector, this->get_objective_multiplier(), multipliers.constraints, result); + const Multipliers& multipliers, double* result, const std::optional& scaling) const { + hessian_model.compute_hessian_vector_product(x, vector, this->get_objective_multiplier(), multipliers.constraints, + result, scaling); } size_t OptimizationProblem::get_number_original_variables() const { diff --git a/uno/optimization/OptimizationProblem.hpp b/uno/optimization/OptimizationProblem.hpp index 4e1c57ce9..89152d866 100644 --- a/uno/optimization/OptimizationProblem.hpp +++ b/uno/optimization/OptimizationProblem.hpp @@ -56,7 +56,7 @@ namespace uno { const Vector& primal_variables, const Multipliers& multipliers, double* hessian_values, const std::optional& scaling) const; virtual void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, - const Multipliers& multipliers, double* result) const; + const Multipliers& multipliers, double* result, const std::optional& scaling) const; [[nodiscard]] size_t get_number_original_variables() const; [[nodiscard]] virtual double variable_lower_bound(size_t variable_index) const;