[WIP] Function scaling

interfaces/AMPL/AMPLModel.cpp

@@ -91,12 +91,13 @@ namespace uno {
 
    double AMPLModel::evaluate_objective(const Vector<double>& x) const {
       fint error_flag = 0;
-      const double result = this->optimization_sense * (*(this->asl)->p.Objval)(this->asl, 0, const_cast<double*>(x.data()), &error_flag);
+      double objective_value = (*(this->asl)->p.Objval)(this->asl, 0, const_cast<double*>(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;
    }
 
    /*

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"

uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp

@@ -69,7 +69,14 @@ 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);
+
+      // 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"));
+      }
+      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);
@@ -157,7 +164,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';
    }
@@ -200,11 +208,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);
 
@@ -227,13 +237,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);
       }

uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp

@@ -5,13 +5,15 @@
 #define UNO_FEASIBILITYRESTORATION_H
 
 #include <memory>
+#include <optional>
 #include "ConstraintRelaxationStrategy.hpp"
 #include "relaxed_problems/l1RelaxedProblem.hpp"
 #include "ingredients/globalization_strategies/GlobalizationStrategy.hpp"
 #include "ingredients/globalization_strategies/l1MeritFunction.hpp"
 #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<InequalityHandlingMethod> feasibility_inequality_handling_method;
       std::unique_ptr<GlobalizationStrategy> optimality_globalization_strategy;
       l1MeritFunction feasibility_globalization_strategy;
+      std::optional<Scaling> 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;

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"
 
@@ -36,10 +37,19 @@ 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, 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);
+
+      // 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"));
+      }
+      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);
    }
@@ -49,7 +59,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();
@@ -68,7 +79,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;
@@ -79,7 +90,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);
    }

uno/ingredients/constraint_relaxation_strategies/NoRelaxation.hpp

@@ -5,11 +5,13 @@
 #define UNO_NORELAXATION_H
 
 #include <memory>
+#include <optional>
 #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<HessianModel> hessian_model;
       std::unique_ptr<InertiaCorrectionStrategy<double>> inertia_correction_strategy;
       l1MeritFunction globalization_strategy;
+      std::optional<Scaling> scaling{};
    };
 } // namespace
 

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>& 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<double>& lagrangian_gradient,
-         const EvaluationSpace& evaluation_space, Iterate& iterate) const {
+         const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional<Scaling>& /*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<double>& primal_variables,
-         const Multipliers& multipliers, double* hessian_values) const {
+         const Multipliers& multipliers, double* hessian_values, const std::optional<Scaling>& 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();
@@ -208,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>& 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.) {

uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp

@@ -24,7 +24,7 @@ namespace uno {
       void evaluate_constraints(Iterate& iterate, Vector<double>& 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>& scaling) const override;
 
       // sparsity patterns of Jacobian and Hessian
       [[nodiscard]] size_t number_jacobian_nonzeros() const override;
@@ -38,11 +38,11 @@ namespace uno {
       // numerical evaluations of Jacobian and Hessian
       void evaluate_constraint_jacobian(Iterate& iterate, double* jacobian_values) const override;
       void evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient,
-         const EvaluationSpace& evaluation_space, Iterate& iterate) const override;
+         const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional<Scaling>& scaling) const override;
       void evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector<double>& primal_variables,
-         const Multipliers& multipliers, double* hessian_values) const override;
+         const Multipliers& multipliers, double* hessian_values, const std::optional<Scaling>& 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>& 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;

uno/ingredients/hessian_models/ExactHessian.cpp

@@ -34,13 +34,20 @@ namespace uno {
    }
 
    void ExactHessian::evaluate_hessian(Statistics& /*statistics*/, const Vector<double>& primal_variables,
-         double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values) {
+         double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values,
+         const std::optional<Scaling>& 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;
    }
 
-   void ExactHessian::compute_hessian_vector_product(const double* x, const double* vector,
-         double objective_multiplier, const Vector<double>& constraint_multipliers, double* result) {
+   void ExactHessian::compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier,
+         const Vector<double>& constraint_multipliers, double* result, const std::optional<Scaling>& 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;
    }

uno/ingredients/hessian_models/ExactHessian.hpp

@@ -22,9 +22,9 @@ namespace uno {
       [[nodiscard]] bool is_positive_definite() const override;
 
       void evaluate_hessian(Statistics& statistics, const Vector<double>& primal_variables, double objective_multiplier,
-         const Vector<double>& constraint_multipliers, double* hessian_values) override;
+         const Vector<double>& constraint_multipliers, double* hessian_values, const std::optional<Scaling>& scaling) override;
       void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier,
-         const Vector<double>& constraint_multipliers, double* result) override;
+         const Vector<double>& constraint_multipliers, double* result, const std::optional<Scaling>& scaling) override;
 
    protected:
       const Model& model;

uno/ingredients/hessian_models/HessianModel.hpp

@@ -5,9 +5,11 @@
 #define UNO_HESSIANMODEL_H
 
 #include <cstddef>
+#include <optional>
 #include <string>
 #include <string_view>
 #include "../interfaces/C/uno_int.h"
+#include "optimization/Scaling.hpp"
 
 namespace uno {
    // forward declarations
@@ -32,9 +34,11 @@ namespace uno {
       [[nodiscard]] virtual bool is_positive_definite() const = 0;
 
       virtual void evaluate_hessian(Statistics& statistics, const Vector<double>& primal_variables,
-         double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values) = 0;
+         double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values,
+         const std::optional<Scaling>& scaling) = 0;
       virtual void compute_hessian_vector_product(const double* x, const double* vector,
-         double objective_multiplier, const Vector<double>& constraint_multipliers, double* result) = 0;
+         double objective_multiplier, const Vector<double>& constraint_multipliers, double* result,
+         const std::optional<Scaling>& scaling) = 0;
    };
 } // namespace
 

uno/ingredients/hessian_models/IdentityHessian.cpp

@@ -38,15 +38,17 @@ namespace uno {
    }
 
    void IdentityHessian::evaluate_hessian(Statistics& /*statistics*/, const Vector<double>& /*primal_variables*/,
-         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* hessian_values) {
+         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* hessian_values,
+         const std::optional<Scaling>& /*scaling*/) {
       DEBUG << "Setting identity Hessian\n";
       for (size_t variable_index: Range(this->number_variables)) {
          hessian_values[variable_index] = 1.;
       }
    }
 
    void IdentityHessian::compute_hessian_vector_product(const double* /*x*/, const double* vector,
-         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* result) {
+         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* result,
+         const std::optional<Scaling>& /*scaling*/) {
       for (size_t variable_index: Range(this->number_variables)) {
          result[variable_index] = vector[variable_index];
       }