Restored Preprocessing::compute_least_square_multipliers (uses a new …

…IdentityMatrix)

uno/ingredients/convexification_strategies/PrimalDualConvexificationStrategy.hpp

@@ -56,7 +56,6 @@ namespace uno {
    void PrimalDualConvexificationStrategy<ElementType>::regularize_matrix(Statistics& statistics,
          DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver, SymmetricMatrix<size_t, ElementType>& matrix,
          size_t size_primal_block, size_t size_dual_block, ElementType dual_regularization_parameter) {
-      DEBUG2 << "Original matrix\n" << matrix << '\n';
       this->primal_regularization = ElementType(0.);
       this->dual_regularization = ElementType(0.);
       DEBUG << "Testing factorization with regularization factors (0, 0)\n";

uno/ingredients/hessian_models/ConvexifiedHessian.hpp

@@ -1,6 +1,9 @@
 // Copyright (c) 2024 Charlie Vanaret
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
+#ifndef UNO_CONVEXIFIEDHESSIAN_H
+#define UNO_CONVEXIFIEDHESSIAN_H
+
 #include <memory>
 #include "HessianModel.hpp"
 
@@ -27,4 +30,6 @@ namespace uno {
 
       void regularize(Statistics& statistics, SymmetricMatrix<size_t, double>& hessian, size_t number_original_variables);
    };
-} // namespace
+} // namespace
+
+#endif // UNO_CONVEXIFIEDHESSIAN_H

uno/ingredients/hessian_models/ExactHessian.hpp

@@ -1,6 +1,9 @@
 // Copyright (c) 2024 Charlie Vanaret
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
+#ifndef UNO_EXACTHESSIAN_H
+#define UNO_EXACTHESSIAN_H
+
 #include "HessianModel.hpp"
 
 namespace uno {
@@ -13,4 +16,6 @@ namespace uno {
       void evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
             const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) override;
    };
-} // namespace
+} // namespace
+
+#endif // UNO_EXACTHESSIAN_H

uno/ingredients/hessian_models/HessianModelFactory.cpp

@@ -6,6 +6,7 @@
 #include "HessianModel.hpp"
 #include "ConvexifiedHessian.hpp"
 #include "ExactHessian.hpp"
+#include "IdentityHessian.hpp"
 #include "ZeroHessian.hpp"
 #include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
 
@@ -23,6 +24,9 @@ namespace uno {
       else if (hessian_model == "zero") {
          return std::make_unique<ZeroHessian>();
       }
+      else if (hessian_model == "identity") {
+         return std::make_unique<IdentityHessian>();
+      }
       throw std::invalid_argument("Hessian model " + hessian_model + " does not exist");
    }
 } // namespace

uno/ingredients/hessian_models/IdentityHessian.cpp

@@ -0,0 +1,24 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#include "IdentityHessian.hpp"
+#include "optimization/OptimizationProblem.hpp"
+#include "linear_algebra/SymmetricMatrix.hpp"
+#include "options/Options.hpp"
+
+namespace uno {
+   // identity Hessian
+   IdentityHessian::IdentityHessian(): HessianModel() {
+   }
+
+   void IdentityHessian::initialize_statistics(Statistics& /*statistics*/, const Options& /*options*/) const { }
+
+   void IdentityHessian::evaluate(Statistics& /*statistics*/, const OptimizationProblem& problem, const Vector<double>& /*primal_variables*/,
+         const Vector<double>& /*constraint_multipliers*/, SymmetricMatrix<size_t, double>& hessian) {
+      // evaluate Lagrangian Hessian
+      hessian.set_dimension(problem.number_variables);
+      for (size_t variable_index: Range(problem.number_variables)) {
+         hessian.insert(1., variable_index, variable_index);
+      }
+   }
+} // namespace

uno/ingredients/hessian_models/IdentityHessian.hpp

@@ -0,0 +1,21 @@
+// Copyright (c) 2024 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#ifndef UNO_IDENTITYHESSIAN_H
+#define UNO_IDENTITYHESSIAN_H
+
+#include "HessianModel.hpp"
+
+namespace uno {
+   // exact Hessian
+   class IdentityHessian : public HessianModel {
+   public:
+      IdentityHessian();
+
+      void initialize_statistics(Statistics& statistics, const Options& options) const override;
+      void evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
+            const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) override;
+   };
+} // namespace
+
+#endif // UNO_IDENTITYHESSIAN_H

uno/ingredients/hessian_models/ZeroHessian.hpp

@@ -1,6 +1,9 @@
 // Copyright (c) 2024 Charlie Vanaret
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
+#ifndef UNO_ZEROHESSIAN_H
+#define UNO_ZEROHESSIAN_H
+
 #include "HessianModel.hpp"
 
 namespace uno {
@@ -13,4 +16,6 @@ namespace uno {
       void evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
             const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) override;
    };
-} // namespace
+} // namespace
+
+#endif // UNO_ZEROHESSIAN_H

uno/ingredients/inequality_handling_methods/inequality_constrained_methods/LPMethod.cpp

@@ -24,7 +24,7 @@ namespace uno {
 
    void LPMethod::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,
          const Multipliers& current_multipliers, Direction& direction, WarmstartInformation& warmstart_information) {
-      LagrangeNewtonSubproblem subproblem(problem, current_iterate, current_multipliers.constraints, *this->hessian_model, this->trust_region_radius);
+      LagrangeNewtonSubproblem subproblem(problem, current_iterate, current_multipliers, *this->hessian_model, this->trust_region_radius);
       this->solver->solve_LP(statistics, subproblem, this->initial_point, direction, warmstart_information);
       InequalityConstrainedMethod::compute_dual_displacements(current_multipliers, direction.multipliers);
       this->number_subproblems_solved++;

uno/ingredients/inequality_handling_methods/inequality_constrained_methods/QPMethod.cpp

@@ -37,7 +37,7 @@ namespace uno {
 
    void QPMethod::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,  const Multipliers& current_multipliers,
          Direction& direction, WarmstartInformation& warmstart_information) {
-      LagrangeNewtonSubproblem subproblem(problem, current_iterate, current_multipliers.constraints, *this->hessian_model, this->trust_region_radius);
+      LagrangeNewtonSubproblem subproblem(problem, current_iterate, current_multipliers, *this->hessian_model, this->trust_region_radius);
       this->solver->solve_QP(statistics, subproblem, this->initial_point, direction, warmstart_information);
       InequalityConstrainedMethod::compute_dual_displacements(current_multipliers, direction.multipliers);
       this->number_subproblems_solved++;

uno/ingredients/inequality_handling_methods/interior_point_methods/PrimalDualInteriorPointMethod.cpp

@@ -41,7 +41,7 @@ namespace uno {
                options.get_double("barrier_push_variable_to_interior_k1"),
                options.get_double("barrier_push_variable_to_interior_k2")
          }),
-         // least_square_multiplier_max_norm(options.get_double("least_square_multiplier_max_norm")),
+         least_square_multiplier_max_norm(options.get_double("least_square_multiplier_max_norm")),
          damping_factor(options.get_double("barrier_damping_factor")),
          l1_constraint_violation_coefficient(options.get_double("l1_constraint_violation_coefficient")) {
    }
@@ -100,7 +100,8 @@ namespace uno {
 
       // compute least-square multipliers
       if (problem.is_constrained()) {
-         this->compute_least_square_multipliers(problem, initial_iterate, initial_iterate.multipliers.constraints);
+         Preprocessing::compute_least_square_multipliers(problem, *this->linear_solver, initial_iterate, initial_iterate.multipliers,
+               this->least_square_multiplier_max_norm);
       }
    }
 
@@ -140,8 +141,7 @@ namespace uno {
 
       // create the barrier reformulation and the subproblem
       PrimalDualInteriorPointProblem barrier_problem(problem, current_multipliers, this->barrier_parameter());
-      LagrangeNewtonSubproblem subproblem(barrier_problem, current_iterate, current_multipliers.constraints, *this->hessian_model,
-            this->trust_region_radius);
+      LagrangeNewtonSubproblem subproblem(barrier_problem, current_iterate, current_multipliers, *this->hessian_model, this->trust_region_radius);
       this->linear_solver->solve_indefinite_system(statistics, subproblem, this->solution, warmstart_information);
 
       assert(direction.status == SubproblemStatus::OPTIMAL && "The primal-dual perturbed subproblem was not solved to optimality");
@@ -211,7 +211,8 @@ namespace uno {
       assert(this->solving_feasibility_problem && "The barrier subproblem did not know it was solving the feasibility problem.");
       this->barrier_parameter_update_strategy.set_barrier_parameter(this->previous_barrier_parameter);
       this->solving_feasibility_problem = false;
-      this->compute_least_square_multipliers(problem, trial_iterate, trial_iterate.multipliers.constraints);
+      Preprocessing::compute_least_square_multipliers(problem, *this->linear_solver, trial_iterate, trial_iterate.multipliers,
+            this->least_square_multiplier_max_norm);
    }
 
    void PrimalDualInteriorPointMethod::set_auxiliary_measure(const Model& model, Iterate& iterate) {
@@ -384,16 +385,6 @@ namespace uno {
       }
    }
 
-   void PrimalDualInteriorPointMethod::compute_least_square_multipliers(const OptimizationProblem& /*problem*/, Iterate& /*iterate*/,
-         Vector<double>& /*constraint_multipliers*/) {
-      /*
-      this->augmented_system.matrix.set_dimension(problem.number_variables + problem.number_constraints);
-      this->augmented_system.matrix.reset();
-      Preprocessing::compute_least_square_multipliers(problem.model, this->augmented_system.matrix, this->augmented_system.rhs, *this->linear_solver,
-            iterate, constraint_multipliers, this->least_square_multiplier_max_norm);
-            */
-   }
-
    void PrimalDualInteriorPointMethod::postprocess_iterate(const OptimizationProblem& problem, Iterate& iterate) {
       // rescale the bound multipliers (Eq. 16 in Ipopt paper)
       for (const size_t variable_index: problem.get_lower_bounded_variables()) {

uno/ingredients/inequality_handling_methods/interior_point_methods/PrimalDualInteriorPointMethod.hpp

@@ -58,7 +58,7 @@ namespace uno {
       double previous_barrier_parameter;
       const double default_multiplier;
       const InteriorPointParameters parameters;
-      // const double least_square_multiplier_max_norm;
+      const double least_square_multiplier_max_norm;
       const double damping_factor; // (Section 3.7 in IPOPT paper)
       const double l1_constraint_violation_coefficient; // (rho in Section 3.3.1 in IPOPT paper)
 
@@ -81,7 +81,6 @@ namespace uno {
             Vector<double>& direction_primals, Multipliers& direction_multipliers);
       void compute_bound_dual_direction(const OptimizationProblem& problem, const Vector<double>& current_primals, const Multipliers& current_multipliers,
             const Vector<double>& primal_direction, Multipliers& direction_multipliers);
-      void compute_least_square_multipliers(const OptimizationProblem& problem, Iterate& iterate, Vector<double>& constraint_multipliers);
    };
 } // namespace
 

uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp

@@ -151,16 +151,18 @@ namespace uno {
       // constraints and Jacobian
       if (warmstart_information.constraints_changed) {
          subproblem.evaluate_constraints(this->constraints);
+      }
+      if (warmstart_information.constraint_jacobian_changed) {
          subproblem.evaluate_constraint_jacobian(this->constraint_jacobian);
       }
 
       // Lagrangian Hessian
-      if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
+      if (warmstart_information.objective_changed || warmstart_information.constraints_changed || warmstart_information.constraint_jacobian_changed) {
          DEBUG << "Evaluating problem Hessian\n";
          subproblem.evaluate_hessian(statistics, this->hessian);
       }
 
-      if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
+      if (warmstart_information.objective_changed || warmstart_information.constraint_jacobian_changed) {
          // form the KKT matrix
          this->augmented_matrix.set_dimension(subproblem.number_variables + subproblem.number_constraints);
          this->augmented_matrix.reset();
@@ -185,7 +187,7 @@ namespace uno {
          this->do_symbolic_analysis(this->augmented_matrix);
          warmstart_information.hessian_sparsity_changed = warmstart_information.jacobian_sparsity_changed = false;
       }
-      if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
+      if (warmstart_information.objective_changed || warmstart_information.constraint_jacobian_changed) {
          DEBUG << "Performing numerical factorization of the augmented matrix\n";
          this->do_numerical_factorization(this->augmented_matrix);
          // regularize

uno/ingredients/subproblems/LagrangeNewtonSubproblem.cpp

@@ -8,7 +8,7 @@
 
 namespace uno {
    LagrangeNewtonSubproblem::LagrangeNewtonSubproblem(const OptimizationProblem& problem, Iterate& current_iterate,
-      const Vector<double>& current_multipliers, HessianModel& hessian_model, double trust_region_radius):
+      const Multipliers& current_multipliers, HessianModel& hessian_model, double trust_region_radius):
             number_variables(problem.number_variables), number_constraints(problem.number_constraints),
             problem(problem), current_iterate(current_iterate), current_multipliers(current_multipliers), hessian_model(hessian_model),
             trust_region_radius(trust_region_radius) { }
@@ -26,7 +26,7 @@ namespace uno {
    }
 
    void LagrangeNewtonSubproblem::evaluate_hessian(Statistics& statistics, SymmetricMatrix<size_t, double>& hessian) {
-      this->hessian_model.evaluate(statistics, this->problem, this->current_iterate.primals, this->current_multipliers, hessian);
+      this->hessian_model.evaluate(statistics, this->problem, this->current_iterate.primals, this->current_multipliers.constraints, hessian);
    }
 
    void LagrangeNewtonSubproblem::compute_lagrangian_gradient(SparseVector<double>& objective_gradient, RectangularMatrix<double>& jacobian,
@@ -40,12 +40,19 @@ namespace uno {
 
       // constraints
       for (size_t constraint_index: Range(this->number_constraints)) {
-         if (this->current_multipliers[constraint_index] != 0.) {
+         if (this->current_multipliers.constraints[constraint_index] != 0.) {
             for (const auto [variable_index, derivative]: jacobian[constraint_index]) {
-               gradient[variable_index] -= this->current_multipliers[constraint_index] * derivative;
+               gradient[variable_index] -= this->current_multipliers.constraints[constraint_index] * derivative;
             }
          }
       }
+
+      /*
+      // bound multipliers
+      for (size_t variable_index: Range(this->number_variables)) {
+         gradient[variable_index] -= (this->current_multipliers.lower_bounds[variable_index] + this->current_multipliers.upper_bounds[variable_index]);
+      }
+       */
    }
 
    double LagrangeNewtonSubproblem::dual_regularization_parameter() const {

uno/ingredients/subproblems/LagrangeNewtonSubproblem.hpp

@@ -13,6 +13,7 @@ namespace uno {
    // forward declarations
    class HessianModel;
    class Iterate;
+   class Multipliers;
    class OptimizationProblem;
    template <typename ElementType>
    class RectangularMatrix;
@@ -29,7 +30,7 @@ namespace uno {
       const size_t number_variables;
       const size_t number_constraints;
 
-      LagrangeNewtonSubproblem(const OptimizationProblem& problem, Iterate& current_iterate, const Vector<double>& current_multipliers,
+      LagrangeNewtonSubproblem(const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers,
             HessianModel& hessian_model, double trust_region_radius);
 
       void evaluate_objective_gradient(SparseVector<double>& objective_gradient);
@@ -48,7 +49,7 @@ namespace uno {
    protected:
       const OptimizationProblem& problem;
       Iterate& current_iterate;
-      const Vector<double>& current_multipliers;
+      const Multipliers& current_multipliers;
       HessianModel& hessian_model;
       // TODO convexification model
       double trust_region_radius;

uno/optimization/WarmstartInformation.cpp

@@ -17,6 +17,7 @@ namespace uno {
    void WarmstartInformation::no_changes() {
       this->objective_changed = false;
       this->constraints_changed = false;
+      this->constraint_jacobian_changed = false;
       this->constraint_bounds_changed = false;
       this->variable_bounds_changed = false;
       this->hessian_sparsity_changed = false;
@@ -26,13 +27,15 @@ namespace uno {
    void WarmstartInformation::iterate_changed() {
       this->objective_changed = true;
       this->constraints_changed = true;
+      this->constraint_jacobian_changed = true;
       this->constraint_bounds_changed = true;
       this->variable_bounds_changed = true;
    }
 
    void WarmstartInformation::whole_problem_changed() {
       this->objective_changed = true;
       this->constraints_changed = true;
+      this->constraint_jacobian_changed = true;
       this->constraint_bounds_changed = true;
       this->variable_bounds_changed = true;
       this->hessian_sparsity_changed = true;
@@ -42,6 +45,7 @@ namespace uno {
    void WarmstartInformation::only_objective_changed() {
       this->objective_changed = true;
       this->constraints_changed = false;
+      this->constraint_jacobian_changed = false;
       this->constraint_bounds_changed = false;
       this->variable_bounds_changed = false;
       this->hessian_sparsity_changed = false;

uno/optimization/WarmstartInformation.hpp

@@ -15,6 +15,7 @@ namespace uno {
 
       bool objective_changed{true};
       bool constraints_changed{true};
+      bool constraint_jacobian_changed{true};
       bool constraint_bounds_changed{true};
       bool variable_bounds_changed{true};
       // bool problem_structure_changed{true};