Merge pull request #104 from cvanaret/improve_factorizations

Reduce the number of symbolic analyses

.github/julia/runtests.jl

@@ -18,11 +18,11 @@ import Uno_jll
 Create a new `AmplNLWriter.Optimizer` object that uses Uno as the backing
 solver.
 """
-function Optimizer(options = String["logger=SILENT"])
+function Optimizer(options = String["logger=SILENT", "unbounded_objective_threshold=-1e15"])
     return AmplNLWriter.Optimizer(Uno_jll.amplexe, options)
 end
 
-Optimizer_LP() = Optimizer(["logger=SILENT", "preset=filterslp", "max_iterations=10000"])
+Optimizer_LP() = Optimizer(["logger=SILENT", "preset=filterslp", "max_iterations=10000", "unbounded_objective_threshold=-1e15"])
 
 # This testset runs https://github.com/jump-dev/MINLPTests.jl
 @testset "MINLPTests" begin

uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp

@@ -3,12 +3,13 @@
 
 #include <cmath>
 #include "PrimalDualInteriorPointSubproblem.hpp"
-#include "optimization/Direction.hpp"
-#include "optimization/Iterate.hpp"
 #include "ingredients/hessian_models/HessianModelFactory.hpp"
 #include "linear_algebra/SparseStorageFactory.hpp"
+#include "linear_algebra/SymmetricIndefiniteLinearSystem.hpp"
 #include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
 #include "solvers/SymmetricIndefiniteLinearSolverFactory.hpp"
+#include "optimization/Direction.hpp"
+#include "optimization/Iterate.hpp"
 #include "optimization/WarmstartInformation.hpp"
 #include "preprocessing/Preprocessing.hpp"
 #include "reformulation/l1RelaxedProblem.hpp"
@@ -123,26 +124,6 @@ namespace uno {
 
    void PrimalDualInteriorPointSubproblem::evaluate_functions(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,
          const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) {
-      // barrier Lagrangian Hessian
-      if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
-         // original Lagrangian Hessian
-         this->hessian_model->evaluate(statistics, problem, current_iterate.primals, current_multipliers.constraints);
-
-         // diagonal barrier terms (grouped by variable)
-         for (size_t variable_index: Range(problem.number_variables)) {
-            double diagonal_barrier_term = 0.;
-            if (is_finite(problem.variable_lower_bound(variable_index))) { // lower bounded
-               const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_lower_bound(variable_index);
-               diagonal_barrier_term += current_multipliers.lower_bounds[variable_index] / distance_to_bound;
-            }
-            if (is_finite(problem.variable_upper_bound(variable_index))) { // upper bounded
-               const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_upper_bound(variable_index);
-               diagonal_barrier_term += current_multipliers.upper_bounds[variable_index] / distance_to_bound;
-            }
-            this->hessian_model->hessian.insert(diagonal_barrier_term, variable_index, variable_index);
-         }
-      }
-
       // barrier objective gradient
       if (warmstart_information.objective_changed) {
          // original objective gradient
@@ -174,6 +155,26 @@ namespace uno {
          problem.evaluate_constraints(current_iterate, this->constraints);
          problem.evaluate_constraint_jacobian(current_iterate, this->constraint_jacobian);
       }
+
+      // barrier Lagrangian Hessian
+      if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
+         // original Lagrangian Hessian
+         this->hessian_model->evaluate(statistics, problem, current_iterate.primals, current_multipliers.constraints);
+
+         // diagonal barrier terms (grouped by variable)
+         for (size_t variable_index: Range(problem.number_variables)) {
+            double diagonal_barrier_term = 0.;
+            if (is_finite(problem.variable_lower_bound(variable_index))) { // lower bounded
+               const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_lower_bound(variable_index);
+               diagonal_barrier_term += current_multipliers.lower_bounds[variable_index] / distance_to_bound;
+            }
+            if (is_finite(problem.variable_upper_bound(variable_index))) { // upper bounded
+               const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_upper_bound(variable_index);
+               diagonal_barrier_term += current_multipliers.upper_bounds[variable_index] / distance_to_bound;
+            }
+            this->hessian_model->hessian.insert(diagonal_barrier_term, variable_index, variable_index);
+         }
+      }
    }
 
    void PrimalDualInteriorPointSubproblem::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,
@@ -199,7 +200,7 @@ namespace uno {
       this->evaluate_functions(statistics, problem, current_iterate, current_multipliers, warmstart_information);
 
       // compute the primal-dual solution
-      this->assemble_augmented_system(statistics, problem, current_multipliers);
+      this->assemble_augmented_system(statistics, problem, current_multipliers, warmstart_information);
       this->augmented_system.solve(*this->linear_solver);
       assert(direction.status == SubproblemStatus::OPTIMAL && "The primal-dual perturbed subproblem was not solved to optimality");
       this->number_subproblems_solved++;
@@ -209,13 +210,14 @@ namespace uno {
    }
 
    void PrimalDualInteriorPointSubproblem::assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem,
-         const Multipliers& current_multipliers) {
+         const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) {
       // assemble, factorize and regularize the augmented matrix
       this->augmented_system.assemble_matrix(this->hessian_model->hessian, this->constraint_jacobian, problem.number_variables, problem.number_constraints);
-      this->augmented_system.factorize_matrix(problem.model, *this->linear_solver);
+      DEBUG << "Testing factorization with regularization factors (0, 0)\n";
+      this->augmented_system.factorize_matrix(*this->linear_solver, warmstart_information);
       const double dual_regularization_parameter = std::pow(this->barrier_parameter(), this->parameters.regularization_exponent);
-      this->augmented_system.regularize_matrix(statistics, problem.model, *this->linear_solver, problem.number_variables, problem.number_constraints,
-            dual_regularization_parameter);
+      this->augmented_system.regularize_matrix(statistics, *this->linear_solver, problem.number_variables, problem.number_constraints,
+            dual_regularization_parameter, warmstart_information);
 
       // check the inertia
       [[maybe_unused]] auto [number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues] = this->linear_solver->get_inertia();

uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.hpp

@@ -79,7 +79,8 @@ namespace uno {
             const Vector<double>& primal_direction, double tau);
       [[nodiscard]] static double dual_fraction_to_boundary(const OptimizationProblem& problem, const Multipliers& current_multipliers,
             Multipliers& direction_multipliers, double tau);
-      void assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem, const Multipliers& current_multipliers);
+      void assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem, const Multipliers& current_multipliers,
+            const WarmstartInformation& warmstart_information);
       void assemble_augmented_rhs(const OptimizationProblem& problem, const Multipliers& current_multipliers);
       void assemble_primal_dual_direction(const OptimizationProblem& problem, const Vector<double>& current_primals, const Multipliers& current_multipliers,
             Vector<double>& direction_primals, Multipliers& direction_multipliers);

uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp

@@ -10,8 +10,9 @@
 #include "RectangularMatrix.hpp"
 #include "ingredients/hessian_models/UnstableRegularization.hpp"
 #include "model/Model.hpp"
-#include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
+#include "optimization/WarmstartInformation.hpp"
 #include "options/Options.hpp"
+#include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
 #include "tools/Statistics.hpp"
 
 namespace uno {
@@ -26,14 +27,13 @@ namespace uno {
             const Options& options);
       void assemble_matrix(const SymmetricMatrix<size_t, double>& hessian, const RectangularMatrix<double>& constraint_jacobian,
             size_t number_variables, size_t number_constraints);
-      void factorize_matrix(const Model& model, DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver);
-      void regularize_matrix(Statistics& statistics, const Model& model, DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver,
-            size_t size_primal_block, size_t size_dual_block, ElementType dual_regularization_parameter);
+      void factorize_matrix(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver, const WarmstartInformation& warmstart_information);
+      void regularize_matrix(Statistics& statistics, DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver,
+            size_t size_primal_block, size_t size_dual_block, ElementType dual_regularization_parameter, const WarmstartInformation& warmstart_information);
       void solve(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver);
       // [[nodiscard]] T get_primal_regularization() const;
 
    protected:
-      size_t number_factorizations{0};
       ElementType primal_regularization{0.};
       ElementType dual_regularization{0.};
       ElementType previous_primal_regularization{0.};
@@ -89,34 +89,36 @@ namespace uno {
    }
 
    template <typename ElementType>
-   void SymmetricIndefiniteLinearSystem<ElementType>::factorize_matrix(const Model& /*model*/,
-         DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver) {
-      if (true || this->number_factorizations == 0) {
+   void SymmetricIndefiniteLinearSystem<ElementType>::factorize_matrix(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver,
+         const WarmstartInformation& warmstart_information) {
+      if (warmstart_information.problem_structure_changed) {
+         DEBUG << "Performing symbolic analysis of the indefinite system\n";
          linear_solver.do_symbolic_analysis(this->matrix);
       }
+      DEBUG << "Performing numerical factorization of the indefinite system\n";
       linear_solver.do_numerical_factorization(this->matrix);
-      this->number_factorizations++;
    }
 
+   // the matrix has been factorized prior to calling this function
    template <typename ElementType>
-   void SymmetricIndefiniteLinearSystem<ElementType>::regularize_matrix(Statistics& statistics, const Model& model,
+   void SymmetricIndefiniteLinearSystem<ElementType>::regularize_matrix(Statistics& statistics,
          DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver, size_t size_primal_block, size_t size_dual_block,
-         ElementType dual_regularization_parameter) {
+         ElementType dual_regularization_parameter, const WarmstartInformation& warmstart_information) {
       DEBUG2 << "Original matrix\n" << this->matrix << '\n';
       this->primal_regularization = ElementType(0.);
       this->dual_regularization = ElementType(0.);
       size_t number_attempts = 1;
-      DEBUG << "Testing factorization with regularization factors (" << this->primal_regularization << ", " << this->dual_regularization << ")\n";
+      DEBUG << "Number of attempts: " << number_attempts << "\n\n";
+
+      auto [number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues] = linear_solver.get_inertia();
+      DEBUG << "Expected inertia  (" << size_primal_block << ", " << size_dual_block << ", 0)\n";
+      DEBUG << "Estimated inertia (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n";
 
-      if (not linear_solver.matrix_is_singular() && linear_solver.number_negative_eigenvalues() == size_dual_block) {
-         DEBUG << "Inertia is good\n";
+      if (number_pos_eigenvalues == size_primal_block && number_neg_eigenvalues == size_dual_block && number_zero_eigenvalues == 0) {
+         DEBUG << "The inertia is correct\n";
          statistics.set("regulariz", this->primal_regularization);
          return;
       }
-      auto [number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues] = linear_solver.get_inertia();
-      DEBUG << "Expected inertia (" << size_primal_block << ", " << size_dual_block << ", 0), ";
-      DEBUG << "got (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n";
-      DEBUG << "Number of attempts: " << number_attempts << "\n\n";
 
       // set the constraint regularization coefficient
       if (linear_solver.matrix_is_singular()) {
@@ -141,19 +143,20 @@ namespace uno {
       while (not good_inertia) {
          DEBUG << "Testing factorization with regularization factors (" << this->primal_regularization << ", " << this->dual_regularization << ")\n";
          DEBUG2 << this->matrix << '\n';
-         this->factorize_matrix(model, linear_solver);
+         this->factorize_matrix(linear_solver, warmstart_information);
          number_attempts++;
+         DEBUG << "Number of attempts: " << number_attempts << "\n";
+
+         std::tie(number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues) = linear_solver.get_inertia();
+         DEBUG << "Expected inertia  (" << size_primal_block << ", " << size_dual_block << ", 0)\n";
+         DEBUG << "Estimated inertia (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n";
 
-         if (not linear_solver.matrix_is_singular() && linear_solver.number_negative_eigenvalues() == size_dual_block) {
+         if (number_pos_eigenvalues == size_primal_block && number_neg_eigenvalues == size_dual_block && number_zero_eigenvalues == 0) {
             good_inertia = true;
-            DEBUG << "Factorization was a success\n";
+            DEBUG << "The inertia is correct\n";
             this->previous_primal_regularization = this->primal_regularization;
          }
          else {
-            std::tie(number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues) = linear_solver.get_inertia();
-            DEBUG << "Expected inertia (" << size_primal_block << ", " << size_dual_block << ", 0), ";
-            DEBUG << "got (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n";
-            DEBUG << "Number of attempts: " << number_attempts << "\n";
             if (this->previous_primal_regularization == 0. || this->threshold_unsuccessful_attempts < number_attempts) {
                this->primal_regularization *= this->primal_regularization_fast_increase_factor;
             }
@@ -188,4 +191,4 @@ namespace uno {
    */
 } // namespace
 
-#endif // UNO_SYMMETRICINDEFINITELINEARSYSTEM_H
+#endif // UNO_SYMMETRICINDEFINITELINEARSYSTEM_H

uno/optimization/WarmstartInformation.cpp

@@ -10,15 +10,15 @@ namespace uno {
       std::cout << "Constraints: " << std::boolalpha << this->constraints_changed << '\n';
       std::cout << "Constraint bounds: " << std::boolalpha << this->constraint_bounds_changed << '\n';
       std::cout << "Variable bounds: " << std::boolalpha << this->variable_bounds_changed << '\n';
-      std::cout << "Problem: " << std::boolalpha << this->problem_changed << '\n';
+      std::cout << "Problem structure: " << std::boolalpha << this->problem_structure_changed << '\n';
    }
 
    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;
+      this->problem_structure_changed = false;
    }
 
    void WarmstartInformation::iterate_changed() {
@@ -33,14 +33,14 @@ namespace uno {
       this->constraints_changed = true;
       this->constraint_bounds_changed = true;
       this->variable_bounds_changed = true;
-      this->problem_changed = true;
+      this->problem_structure_changed = true;
    }
 
    void WarmstartInformation::only_objective_changed() {
       this->objective_changed = true;
       this->constraints_changed = false;
       this->constraint_bounds_changed = false;
       this->variable_bounds_changed = false;
-      this->problem_changed = false;
+      this->problem_structure_changed = false;
    }
 } // namespace

uno/optimization/WarmstartInformation.hpp

@@ -10,7 +10,7 @@ namespace uno {
       bool constraints_changed{true};
       bool constraint_bounds_changed{true};
       bool variable_bounds_changed{true};
-      bool problem_changed{true};
+      bool problem_structure_changed{true};
 
       void display() const;
       void no_changes();

uno/solvers/BQPD/BQPDSolver.cpp

@@ -165,8 +165,8 @@ namespace uno {
 
    BQPDMode BQPDSolver::determine_mode(const WarmstartInformation& warmstart_information) const {
       BQPDMode mode = (this->number_calls == 0) ? BQPDMode::ACTIVE_SET_EQUALITIES : BQPDMode::USER_DEFINED;
-      // if problem changed, use cold start
-      if (warmstart_information.problem_changed) {
+      // if problem structure changed, use cold start
+      if (warmstart_information.problem_structure_changed) {
          mode = BQPDMode::ACTIVE_SET_EQUALITIES;
       }
       // if only the variable bounds changed, reuse the active set estimate and the Jacobian information

uno/solvers/MA57/MA57Solver.cpp

@@ -165,7 +165,7 @@ namespace uno {
       // build the internal matrix representation
       this->row_indices.clear();
       this->column_indices.clear();
-      for (const auto [row_index, column_index, element]: matrix) {
+      for (const auto [row_index, column_index, _]: matrix) {
          this->row_indices.emplace_back(static_cast<int>(row_index + this->fortran_shift));
          this->column_indices.emplace_back(static_cast<int>(column_index + this->fortran_shift));
       }