[WIP] l1 relaxation (based on SQuID method)

uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategyFactory.cpp

@@ -4,6 +4,7 @@
 #include <string>
 #include "ConstraintRelaxationStrategyFactory.hpp"
 #include "FeasibilityRestoration.hpp"
+#include "l1Relaxation.hpp"
 #include "UnconstrainedStrategy.hpp"
 #include "options/Options.hpp"
 #include "tools/Logger.hpp"
@@ -20,6 +21,9 @@ namespace uno {
       if (constraint_relaxation_type == "feasibility_restoration") {
          return std::make_unique<FeasibilityRestoration>(options);
       }
+      else if (constraint_relaxation_type == "l1_relaxation") {
+         return std::make_unique<l1Relaxation>(options);
+      }
       throw std::invalid_argument("ConstraintRelaxationStrategy " + constraint_relaxation_type + " is not supported");
    }
 } // namespace

uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp

@@ -0,0 +1,75 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#ifndef UNO_L1RELAXATION_H
+#define UNO_L1RELAXATION_H
+
+#include <memory>
+#include "ConstraintRelaxationStrategy.hpp"
+#include "relaxed_problems/l1RelaxedProblem.hpp"
+#include "ingredients/hessian_models/HessianModel.hpp"
+#include "ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp"
+
+namespace uno {
+   struct l1RelaxationParameters {
+      double beta;
+      double theta;
+      double kappa_rho;
+      double kappa_lambda;
+      double epsilon;
+      double omega;
+      double delta;
+   };
+
+   class l1Relaxation : public ConstraintRelaxationStrategy {
+   public:
+      explicit l1Relaxation(const Options& options);
+      ~l1Relaxation() override = default;
+
+      void initialize(Statistics& statistics, const Model& model, Iterate& initial_iterate, Direction& direction,
+         double trust_region_radius, const Options& options) override;
+
+      // direction computation
+      void compute_feasible_direction(Statistics& statistics, GlobalizationStrategy& globalization_strategy, Iterate& current_iterate,
+         Direction& direction, double trust_region_radius, WarmstartInformation& warmstart_information) override;
+      [[nodiscard]] bool solving_feasibility_problem() const override;
+      void switch_to_feasibility_problem(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
+         Iterate& current_iterate, double trust_region_radius, WarmstartInformation& warmstart_information) override;
+
+      // trial iterate acceptance
+      [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
+         double trust_region_radius, const Model& model, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,
+         double step_length, WarmstartInformation& warmstart_information, UserCallbacks& user_callbacks) override;
+      [[nodiscard]] SolutionStatus check_termination(const Model& model, Iterate& iterate) override;
+
+      [[nodiscard]] std::string get_name() const override;
+      [[nodiscard]] size_t get_number_subproblems_solved() const override;
+
+   protected:
+      std::unique_ptr<l1RelaxedProblem> relaxed_problem{};
+      std::unique_ptr<InequalityHandlingMethod> inequality_handling_method;
+      std::unique_ptr<HessianModel> hessian_model{};
+      std::unique_ptr<InertiaCorrectionStrategy<double>> inertia_correction_strategy;
+      // feasibility problem
+      std::unique_ptr<const l1RelaxedProblem> feasibility_problem{};
+      std::unique_ptr<InequalityHandlingMethod> feasibility_inequality_handling_method;
+      std::unique_ptr<HessianModel> feasibility_hessian_model{};
+      Multipliers feasibility_multipliers;
+      double penalty_parameter;
+      const double tolerance;
+      const l1RelaxationParameters parameters;
+
+      // auxiliary functions with the notations of the paper
+      [[nodiscard]] double v(const Iterate& current_iterate) const;
+      [[nodiscard]] double l(const Direction& direction, const Iterate& current_iterate) const;
+      [[nodiscard]] double delta_l(const Direction& direction, const Iterate& current_iterate) const;
+      [[nodiscard]] double Ropt(Iterate& current_iterate, double objective_multiplier, const Multipliers& multipliers) const;
+      [[nodiscard]] double Rinf(Iterate& current_iterate, const Multipliers& multipliers) const;
+      [[nodiscard]] double delta_m(const Direction& direction, const Iterate& current_iterate, double objective_multiplier) const;
+      [[nodiscard]] double compute_zeta(const Direction& direction, const Iterate& current_iterate) const;
+      [[nodiscard]] double compute_w(const Direction& feasibility_direction, const Direction& optimality_direction, const Iterate& current_iterate);
+      void check_exact_relaxation(const Iterate& iterate) const;
+   };
+} // namespace
+
+#endif //UNO_L1RELAXATION_H

uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp

@@ -27,16 +27,20 @@ namespace uno {
          constraint_violation_coefficient(constraint_violation_coefficient) {
    }
 
+   void l1RelaxedProblem::set_objective_multiplier(double new_objective_multiplier) {
+      this->objective_multiplier = new_objective_multiplier;
+   }
+
    double l1RelaxedProblem::get_objective_multiplier() const {
       return this->objective_multiplier;
    }
 
-   void l1RelaxedProblem::set_proximal_coefficient(double proximal_coefficient) {
-      this->proximal_coefficient = proximal_coefficient;
+   void l1RelaxedProblem::set_proximal_coefficient(double new_proximal_coefficient) {
+      this->proximal_coefficient = new_proximal_coefficient;
    }
 
-   void l1RelaxedProblem::set_proximal_center(double* proximal_center) {
-      this->proximal_center = proximal_center;
+   void l1RelaxedProblem::set_proximal_center(double* new_proximal_center) {
+      this->proximal_center = new_proximal_center;
    }
 
    void l1RelaxedProblem::evaluate_constraints(Iterate& iterate, Vector<double>& constraints) const {

uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp

@@ -16,9 +16,10 @@ namespace uno {
          bool relax_linear_constraints);
       ~l1RelaxedProblem() override = default;
 
+      void set_objective_multiplier(double new_objective_multiplier);
       [[nodiscard]] double get_objective_multiplier() const override;
-      void set_proximal_coefficient(double proximal_coefficient);
-      void set_proximal_center(double* proximal_center);
+      void set_proximal_coefficient(double new_proximal_coefficient);
+      void set_proximal_center(double* new_proximal_center);
 
       // constraint evaluations
       void evaluate_constraints(Iterate& iterate, Vector<double>& constraints) const override;
@@ -56,7 +57,7 @@ namespace uno {
       [[nodiscard]] const Collection<size_t>& get_dual_regularization_constraints() const override;
 
       [[nodiscard]] SolutionStatus check_first_order_convergence(const Iterate& current_iterate, double primal_tolerance,
-         double dual_tolerance) const;
+         double dual_tolerance) const override;
 
       void set_elastic_variable_values(Iterate& iterate, const std::function<void(Iterate&, size_t, size_t, double)>& elastic_setting_function) const;
 
@@ -68,7 +69,7 @@ namespace uno {
    protected:
       ElasticVariables elastic_variables;
       const size_t number_elastic_variables;
-      const double objective_multiplier;
+      double objective_multiplier;
       const double constraint_violation_coefficient;
       double proximal_coefficient{INF<double>};
       double* proximal_center{};

uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp

@@ -49,7 +49,6 @@ namespace uno {
       virtual void evaluate_constraint_jacobian(Iterate& iterate) = 0;
       virtual void compute_constraint_jacobian_vector_product(const Vector<double>& vector, Vector<double>& result) const = 0;
       virtual void compute_constraint_jacobian_transposed_vector_product(const Vector<double>& vector, Vector<double>& result) const = 0;
-      [[nodiscard]] virtual double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector<double>& vector) const = 0;
 
       // progress measures
       [[nodiscard]] virtual bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,

uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp

@@ -100,11 +100,6 @@ namespace uno {
       evaluation_space.compute_constraint_jacobian_transposed_vector_product(vector, result);
    }
 
-   double InequalityConstrainedMethod::compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector<double>& vector) const {
-      const auto& evaluation_space = this->solver->get_evaluation_space();
-      return evaluation_space.compute_hessian_quadratic_product(subproblem, vector);
-   }
-
    // compute dual *displacements*
    // because of the way we form LPs/QPs, we get the new *multipliers* back from the solver. To get the dual displacements/direction,
    // we need to subtract the current multipliers

uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp

@@ -35,7 +35,6 @@ namespace uno {
       void evaluate_constraint_jacobian(Iterate& iterate) override;
       void compute_constraint_jacobian_vector_product(const Vector<double>& vector, Vector<double>& result) const override;
       void compute_constraint_jacobian_transposed_vector_product(const Vector<double>& vector, Vector<double>& result) const override;
-      [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector<double>& vector) const override;
 
       // acceptance
       [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,

uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp

@@ -41,7 +41,6 @@ namespace uno {
       void evaluate_constraint_jacobian(Iterate& iterate) override;
       void compute_constraint_jacobian_vector_product(const Vector<double>& vector, Vector<double>& result) const override;
       void compute_constraint_jacobian_transposed_vector_product(const Vector<double>& vector, Vector<double>& result) const override;
-      [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector<double>& vector) const override;
 
       // acceptance
       [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
@@ -259,12 +258,6 @@ namespace uno {
       evaluation_space.compute_constraint_jacobian_transposed_vector_product(vector, result);
    }
 
-   template <typename BarrierProblem>
-   double InteriorPointMethod<BarrierProblem>::compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector<double>& vector) const {
-      const auto& evaluation_space = this->linear_solver->get_evaluation_space();
-      return evaluation_space.compute_hessian_quadratic_product(subproblem, vector);
-   }
-
    template <typename BarrierProblem>
    bool InteriorPointMethod<BarrierProblem>::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
          HessianModel& hessian_model, InertiaCorrectionStrategy<double>& inertia_correction_strategy, double trust_region_radius,

uno/optimization/OptimizationProblem.hpp

@@ -79,7 +79,7 @@ namespace uno {
       [[nodiscard]] virtual double complementarity_error(const Vector<double>& primals, const Vector<double>& constraints,
          const Multipliers& multipliers, double shift_value, Norm residual_norm) const;
 
-      [[nodiscard]] SolutionStatus check_first_order_convergence(const Iterate& current_iterate, double primal_tolerance,
+      [[nodiscard]] virtual SolutionStatus check_first_order_convergence(const Iterate& current_iterate, double primal_tolerance,
          double dual_tolerance) const;
 
       // progress measures

uno/options/Presets.cpp

@@ -133,6 +133,19 @@ namespace uno {
          preset_options.set_bool("switch_to_optimality_requires_linearized_feasibility", true);
          preset_options.set_bool("protect_actual_reduction_against_roundoff", false);
       }
+      else if (preset_name == "squid") {
+         preset_options.set_string("constraint_relaxation_strategy", "l1_relaxation");
+         preset_options.set_string("inequality_handling_method", "inequality_constrained");
+         preset_options.set_string("hessian_model", "exact");
+         preset_options.set_string("inertia_correction_strategy", "primal");
+         preset_options.set_string("globalization_mechanism", "LS");
+         preset_options.set_string("globalization_strategy", "l1_merit");
+         preset_options.set_string("progress_norm", "L1");
+         preset_options.set_string("residual_norm", "INF");
+         preset_options.set_double("l1_constraint_violation_coefficient", 1.);
+         preset_options.set_double("primal_tolerance", 1e-6);
+         preset_options.set_double("dual_tolerance", 1e-6);
+      }
       else {
          throw std::runtime_error("The preset " + preset_name + " is not known");
       }