Merge pull request #83 from cvanaret/TR_duals_fix

TR strategy: properly reset the bounds duals and set the warm-start information

uno/Uno.cpp

@@ -11,6 +11,7 @@
 #include "linear_algebra/Vector.hpp"
 #include "model/Model.hpp"
 #include "optimization/Iterate.hpp"
+#include "optimization/WarmstartInformation.hpp"
 #include "solvers/QPSolverFactory.hpp"
 #include "solvers/LPSolverFactory.hpp"
 #include "solvers/SymmetricIndefiniteLinearSolverFactory.hpp"
@@ -32,11 +33,12 @@ namespace uno {
    void Uno::solve(const Model& model, Iterate& current_iterate, const Options& options) {
       Timer timer{};
       Statistics statistics = Uno::create_statistics(model, options);
+      WarmstartInformation warmstart_information{};
+      warmstart_information.whole_problem_changed();
 
       try {
          // use the initial primal-dual point to initialize the strategies and generate the initial iterate
          this->initialize(statistics, current_iterate, options);
-         current_iterate.status = TerminationStatus::NOT_OPTIMAL;
          // allocate the trial iterate once and for all here
          Iterate trial_iterate(current_iterate);
 
@@ -51,7 +53,8 @@ namespace uno {
                DEBUG << "### Outer iteration " << major_iterations << '\n';
 
                // compute an acceptable iterate by solving a subproblem at the current point
-               this->globalization_mechanism.compute_next_iterate(statistics, model, current_iterate, trial_iterate);
+               warmstart_information.iterate_changed();
+               this->globalization_mechanism.compute_next_iterate(statistics, model, current_iterate, trial_iterate, warmstart_information);
                termination = this->termination_criteria(trial_iterate.status, major_iterations, timer.get_duration());
                // the trial iterate becomes the current iterate for the next iteration
                std::swap(current_iterate, trial_iterate);
@@ -81,6 +84,7 @@ namespace uno {
       this->globalization_mechanism.initialize(statistics, current_iterate, options);
       options.print_used();
       if (Logger::level == INFO) statistics.print_current_line();
+      current_iterate.status = TerminationStatus::NOT_OPTIMAL;
    }
 
    Statistics Uno::create_statistics(const Model& model, const Options& options) {

uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.hpp

@@ -46,7 +46,7 @@ namespace uno {
       void compute_feasible_direction(Statistics& statistics, Iterate& current_iterate, Direction& direction,
             const Vector<double>& initial_point, WarmstartInformation& warmstart_information);
       [[nodiscard]] virtual bool solving_feasibility_problem() const = 0;
-      virtual void switch_to_feasibility_problem(Statistics& statistics, Iterate& current_iterate) = 0;
+      virtual void switch_to_feasibility_problem(Statistics& statistics, Iterate& current_iterate, WarmstartInformation& warmstart_information) = 0;
 
       // trial iterate acceptance
       [[nodiscard]] virtual bool is_iterate_acceptable(Statistics& statistics, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,

uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp

@@ -71,17 +71,16 @@ namespace uno {
                // switch to the feasibility problem, starting from the current direction
                statistics.set("status", "infeas. subproblem");
                DEBUG << "/!\\ The subproblem is infeasible\n";
-               this->switch_to_feasibility_problem(statistics, current_iterate);
-               warmstart_information.set_cold_start();
+               this->switch_to_feasibility_problem(statistics, current_iterate, warmstart_information);
                this->subproblem->set_initial_point(direction.primals);
             }
             else {
+               warmstart_information.no_changes();
                return;
             }
          }
          catch (const UnstableRegularization&) {
-            this->switch_to_feasibility_problem(statistics, current_iterate);
-            warmstart_information.set_cold_start();
+            this->switch_to_feasibility_problem(statistics, current_iterate, warmstart_information);
          }
       }
 
@@ -99,7 +98,8 @@ namespace uno {
    }
 
    // precondition: this->current_phase == Phase::OPTIMALITY
-   void FeasibilityRestoration::switch_to_feasibility_problem(Statistics& statistics, Iterate& current_iterate) {
+   void FeasibilityRestoration::switch_to_feasibility_problem(Statistics& statistics, Iterate& current_iterate,
+         WarmstartInformation& warmstart_information) {
       DEBUG << "Switching from optimality to restoration phase\n";
       this->current_phase = Phase::FEASIBILITY_RESTORATION;
       this->globalization_strategy->notify_switch_to_feasibility(current_iterate.progress);
@@ -114,6 +114,7 @@ namespace uno {
       DEBUG2 << "Current iterate:\n" << current_iterate << '\n';
 
       if (Logger::level == INFO) statistics.print_current_line();
+      warmstart_information.whole_problem_changed();
    }
 
    void FeasibilityRestoration::solve_subproblem(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,
@@ -142,7 +143,7 @@ namespace uno {
 
       this->subproblem->exit_feasibility_problem(this->optimality_problem, trial_iterate);
       // set a cold start in the subproblem solver
-      warmstart_information.set_cold_start();
+      warmstart_information.whole_problem_changed();
    }
 
    bool FeasibilityRestoration::is_iterate_acceptable(Statistics& statistics, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,
@@ -156,6 +157,9 @@ namespace uno {
       if (this->current_phase == Phase::FEASIBILITY_RESTORATION && this->can_switch_to_optimality_phase(current_iterate, trial_iterate, direction, step_length)) {
          this->switch_to_optimality_phase(current_iterate, trial_iterate, warmstart_information);
       }
+      else {
+         warmstart_information.no_changes();
+      }
 
       bool accept_iterate = false;
       if (direction.norm == 0.) {
@@ -170,9 +174,6 @@ namespace uno {
          accept_iterate = this->globalization_strategy->is_iterate_acceptable(statistics, current_iterate.progress, trial_iterate.progress,
                predicted_reduction, this->current_problem().get_objective_multiplier());
       }
-      if (accept_iterate) {
-         // this->set_dual_residuals_statistics(statistics, trial_iterate);
-      }
       ConstraintRelaxationStrategy::set_progress_statistics(statistics, trial_iterate);
       return accept_iterate;
    }

uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp

@@ -25,7 +25,7 @@ namespace uno {
       // direction computation
       void compute_feasible_direction(Statistics& statistics, Iterate& current_iterate, Direction& direction, WarmstartInformation& warmstart_information) override;
       [[nodiscard]] bool solving_feasibility_problem() const override;
-      void switch_to_feasibility_problem(Statistics& statistics, Iterate& current_iterate) override;
+      void switch_to_feasibility_problem(Statistics& statistics, Iterate& current_iterate, WarmstartInformation& warmstart_information) override;
 
       // trial iterate acceptance
       [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,

uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp

@@ -79,7 +79,8 @@ namespace uno {
       return (this->penalty_parameter == 0.);
    }
 
-   void l1Relaxation::switch_to_feasibility_problem(Statistics& /*statistics*/, Iterate& /*current_iterate*/) {
+   void l1Relaxation::switch_to_feasibility_problem(Statistics& /*statistics*/, Iterate& /*current_iterate*/,
+         WarmstartInformation& /*warmstart_information*/) {
       throw std::runtime_error("l1Relaxation::switch_to_feasibility_problem is not implemented");
    }
 

uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp

@@ -32,7 +32,7 @@ namespace uno {
       void compute_feasible_direction(Statistics& statistics, Iterate& current_iterate, Direction& direction,
             WarmstartInformation& warmstart_information) override;
       [[nodiscard]] bool solving_feasibility_problem() const override;
-      void switch_to_feasibility_problem(Statistics& statistics, Iterate& current_iterate) override;
+      void switch_to_feasibility_problem(Statistics& statistics, Iterate& current_iterate, WarmstartInformation& warmstart_information) override;
 
       // trial iterate acceptance
       [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,

uno/ingredients/globalization_mechanisms/BacktrackingLineSearch.cpp

@@ -30,9 +30,8 @@ namespace uno {
       this->constraint_relaxation_strategy.initialize(statistics, initial_iterate, options);
    }
 
-   void BacktrackingLineSearch::compute_next_iterate(Statistics& statistics, const Model& model, Iterate& current_iterate, Iterate& trial_iterate) {
-      WarmstartInformation warmstart_information{};
-      warmstart_information.set_hot_start();
+   void BacktrackingLineSearch::compute_next_iterate(Statistics& statistics, const Model& model, Iterate& current_iterate, Iterate& trial_iterate,
+         WarmstartInformation& warmstart_information) {
       DEBUG2 << "Current iterate\n" << current_iterate << '\n';
 
       this->constraint_relaxation_strategy.compute_feasible_direction(statistics, current_iterate, this->direction, warmstart_information);
@@ -59,6 +58,9 @@ namespace uno {
                   // scale or not the constraint dual direction with the LS step length
                   this->scale_duals_with_step_length ? step_length : 1.);
 
+            // let the constraint relaxation strategy determine which quantities change
+            warmstart_information.no_changes();
+
             is_acceptable = this->constraint_relaxation_strategy.is_iterate_acceptable(statistics, current_iterate, trial_iterate, this->direction,
                   step_length, warmstart_information);
             this->set_statistics(statistics, trial_iterate, this->direction, step_length, number_iterations);
@@ -89,8 +91,7 @@ namespace uno {
                }
                // switch to solving the feasibility problem
                statistics.set("status", "small LS step length");
-               this->constraint_relaxation_strategy.switch_to_feasibility_problem(statistics, current_iterate);
-               warmstart_information.set_cold_start();
+               this->constraint_relaxation_strategy.switch_to_feasibility_problem(statistics, current_iterate, warmstart_information);
                this->constraint_relaxation_strategy.compute_feasible_direction(statistics, current_iterate, this->direction, this->direction.primals,
                      warmstart_information);
                BacktrackingLineSearch::check_unboundedness(this->direction);

uno/ingredients/globalization_mechanisms/BacktrackingLineSearch.hpp

@@ -15,7 +15,8 @@ namespace uno {
       BacktrackingLineSearch(ConstraintRelaxationStrategy& constraint_relaxation_strategy, const Options& options);
 
       void initialize(Statistics& statistics, Iterate& initial_iterate, const Options& options) override;
-      void compute_next_iterate(Statistics& statistics, const Model& model, Iterate& current_iterate, Iterate& trial_iterate) override;
+      void compute_next_iterate(Statistics& statistics, const Model& model, Iterate& current_iterate, Iterate& trial_iterate,
+            WarmstartInformation& warmstart_information) override;
 
    private:
       const double backtracking_ratio;

uno/ingredients/globalization_mechanisms/GlobalizationMechanism.hpp

@@ -13,14 +13,16 @@ namespace uno {
    class Model;
    class Options;
    class Statistics;
+   struct WarmstartInformation;
 
    class GlobalizationMechanism {
    public:
       explicit GlobalizationMechanism(ConstraintRelaxationStrategy& constraint_relaxation_strategy);
       virtual ~GlobalizationMechanism() = default;
 
       virtual void initialize(Statistics& statistics, Iterate& initial_iterate, const Options& options) = 0;
-      virtual void compute_next_iterate(Statistics& statistics, const Model& model, Iterate& current_iterate, Iterate& trial_iterate) = 0;
+      virtual void compute_next_iterate(Statistics& statistics, const Model& model, Iterate& current_iterate, Iterate& trial_iterate,
+            WarmstartInformation& warmstart_information) = 0;
 
       [[nodiscard]] size_t get_hessian_evaluation_count() const;
       [[nodiscard]] size_t get_number_subproblems_solved() const;

uno/ingredients/globalization_mechanisms/TrustRegionStrategy.cpp

@@ -38,9 +38,8 @@ namespace uno {
       this->constraint_relaxation_strategy.initialize(statistics, initial_iterate, options);
    }
 
-   void TrustRegionStrategy::compute_next_iterate(Statistics& statistics, const Model& model, Iterate& current_iterate, Iterate& trial_iterate) {
-      WarmstartInformation warmstart_information{};
-      warmstart_information.set_hot_start();
+   void TrustRegionStrategy::compute_next_iterate(Statistics& statistics, const Model& model, Iterate& current_iterate, Iterate& trial_iterate,
+         WarmstartInformation& warmstart_information) {
       DEBUG2 << "Current iterate\n" << current_iterate << '\n';
 
       size_t number_iterations = 0;
@@ -64,20 +63,23 @@ namespace uno {
                statistics.set("status", "unbounded subproblem");
                if (Logger::level == INFO) statistics.print_current_line();
                this->decrease_radius_aggressively();
-               warmstart_information.set_cold_start();
+               warmstart_information.whole_problem_changed();
             }
             else if (this->direction.status == SubproblemStatus::ERROR) {
                this->set_statistics(statistics, this->direction);
                statistics.set("status", "solver error");
                if (Logger::level == INFO) statistics.print_current_line();
                this->decrease_radius();
-               warmstart_information.set_cold_start();
+               warmstart_information.whole_problem_changed();
             }
             else {
                // take full primal-dual step
                GlobalizationMechanism::assemble_trial_iterate(model, current_iterate, trial_iterate, this->direction, 1., 1.);
                this->reset_active_trust_region_multipliers(model, this->direction, trial_iterate);
 
+               // let the constraint relaxation strategy and the radius update rule determine which quantities change
+               warmstart_information.no_changes();
+
                is_acceptable = this->is_iterate_acceptable(statistics, current_iterate, trial_iterate, this->direction, warmstart_information);
                if (is_acceptable) {
                   this->constraint_relaxation_strategy.set_dual_residuals_statistics(statistics, trial_iterate);
@@ -86,8 +88,7 @@ namespace uno {
                }
                else {
                   this->decrease_radius(this->direction.norm);
-                  // after the first iteration, only the variable bounds are updated
-                  warmstart_information.only_variable_bounds_changed();
+                  warmstart_information.variable_bounds_changed = true;
                }
                if (Logger::level == INFO) statistics.print_current_line();
             }
@@ -97,7 +98,7 @@ namespace uno {
             statistics.set("status", "eval. error");
             if (Logger::level == INFO) statistics.print_current_line();
             this->decrease_radius();
-            warmstart_information.set_cold_start();
+            warmstart_information.whole_problem_changed();
          }
          if (not is_acceptable && this->radius < this->minimum_radius) {
             throw std::runtime_error("Small trust-region radius");
@@ -107,25 +108,24 @@ namespace uno {
 
    void TrustRegionStrategy::reset_active_trust_region_multipliers(const Model& model, const Direction& direction, Iterate& trial_iterate) const {
       assert(0 < this->radius && "The trust-region radius should be positive");
-      // set multipliers for bound constraints active at trust region to 0 (except if one of the original bounds is active)
-      for (size_t variable_index: direction.active_bounds.at_lower_bound) {
-         if (variable_index < model.number_variables && std::abs(direction.primals[variable_index] + this->radius) <= this->activity_tolerance &&
-             this->activity_tolerance < std::abs(trial_iterate.primals[variable_index] - model.variable_lower_bound(variable_index))) {
+      // reset multipliers for bound constraints active at trust region (except if one of the original bounds is active)
+      for (size_t variable_index: Range(model.number_variables)) {
+         if (std::abs(direction.primals[variable_index] + this->radius) <= this->activity_tolerance &&
+               this->activity_tolerance < std::abs(trial_iterate.primals[variable_index] - model.variable_lower_bound(variable_index))) {
             trial_iterate.multipliers.lower_bounds[variable_index] = 0.;
+            trial_iterate.feasibility_multipliers.lower_bounds[variable_index] = 0.;
          }
-      }
-      for (size_t variable_index: direction.active_bounds.at_upper_bound) {
-         if (variable_index < model.number_variables && std::abs(direction.primals[variable_index] - this->radius) <= this->activity_tolerance &&
-             this->activity_tolerance < std::abs(model.variable_upper_bound(variable_index) - trial_iterate.primals[variable_index])) {
+         if (std::abs(direction.primals[variable_index] - this->radius) <= this->activity_tolerance &&
+               this->activity_tolerance < std::abs(model.variable_upper_bound(variable_index) - trial_iterate.primals[variable_index])) {
             trial_iterate.multipliers.upper_bounds[variable_index] = 0.;
+            trial_iterate.feasibility_multipliers.upper_bounds[variable_index] = 0.;
          }
       }
    }
 
    // the trial iterate is accepted by the constraint relaxation strategy or if the step is small and we cannot switch to solving the feasibility problem
    bool TrustRegionStrategy::is_iterate_acceptable(Statistics& statistics, Iterate& current_iterate, Iterate& trial_iterate,
          const Direction& direction, WarmstartInformation& warmstart_information) {
-      // direction.primal_dual_step_length is usually 1, can be lower if reduced by fraction-to-boundary rule
       bool accept_iterate = this->constraint_relaxation_strategy.is_iterate_acceptable(statistics, current_iterate, trial_iterate, direction, 1.,
             warmstart_information);
       this->set_statistics(statistics, trial_iterate, direction);

uno/ingredients/globalization_mechanisms/TrustRegionStrategy.hpp

@@ -12,7 +12,8 @@ namespace uno {
       TrustRegionStrategy(ConstraintRelaxationStrategy& constraint_relaxation_strategy, const Options& options);
 
       void initialize(Statistics& statistics, Iterate& initial_iterate, const Options& options) override;
-      void compute_next_iterate(Statistics& statistics, const Model& model, Iterate& current_iterate, Iterate& trial_iterate) override;
+      void compute_next_iterate(Statistics& statistics, const Model& model, Iterate& current_iterate, Iterate& trial_iterate,
+            WarmstartInformation& warmstart_information) override;
 
    private:
       double radius; /*!< Current trust region radius */

uno/optimization/WarmstartInformation.cpp

@@ -13,20 +13,28 @@ namespace uno {
       std::cout << "Problem: " << std::boolalpha << this->problem_changed << '\n';
    }
 
-   void WarmstartInformation::set_cold_start() {
+   void WarmstartInformation::no_changes() {
+      this->objective_changed = false;
+      this->constraints_changed = false;
+      this->constraint_bounds_changed = false;
+      this->variable_bounds_changed = false;
+      this->problem_changed = false;
+   }
+
+   void WarmstartInformation::iterate_changed() {
       this->objective_changed = true;
       this->constraints_changed = true;
       this->constraint_bounds_changed = true;
       this->variable_bounds_changed = true;
-      this->problem_changed = true;
+      this->problem_changed = false;
    }
 
-   void WarmstartInformation::set_hot_start() {
+   void WarmstartInformation::whole_problem_changed() {
       this->objective_changed = true;
       this->constraints_changed = true;
       this->constraint_bounds_changed = true;
       this->variable_bounds_changed = true;
-      this->problem_changed = false;
+      this->problem_changed = true;
    }
 
    void WarmstartInformation::only_objective_changed() {
@@ -36,12 +44,4 @@ namespace uno {
       this->variable_bounds_changed = false;
       this->problem_changed = false;
    }
-
-   void WarmstartInformation::only_variable_bounds_changed() {
-      this->objective_changed = false;
-      this->constraints_changed = false;
-      this->constraint_bounds_changed = false;
-      this->variable_bounds_changed = true;
-      this->problem_changed = false;
-   }
 } // namespace