[WIP] Make BQPD matrix-free

bindings/AMPL/AMPLModel.cpp

@@ -101,8 +101,8 @@ namespace uno {
       fint error_flag = 0;
       // prevent ASL to crash by catching all evaluation errors
       Jmp_buf err_jmp_uno;
-      asl->i.err_jmp_ = &err_jmp_uno;
-      asl->i.err_jmp1_ = &err_jmp_uno;
+      this->asl->i.err_jmp_ = &err_jmp_uno;
+      this->asl->i.err_jmp1_ = &err_jmp_uno;
       if (setjmp(err_jmp_uno.jb)) {
          error_flag = 1;
       }
@@ -207,6 +207,31 @@ namespace uno {
       //this->asl->i.x_known = 0;
    }
 
+   void AMPLModel::compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+         Vector<double>& result) const {
+      fint error_flag = 0;
+      // prevent ASL to crash by catching all evaluation errors
+      Jmp_buf err_jmp_uno;
+      this->asl->i.err_jmp_ = &err_jmp_uno;
+      this->asl->i.err_jmp1_ = &err_jmp_uno;
+      if (setjmp(err_jmp_uno.jb)) {
+         error_flag = 1;
+      }
+
+      // scale by the objective sign
+      objective_multiplier *= this->objective_sign;
+      const int objective_number = -1;
+      // flip the signs of the multipliers: in AMPL, the Lagrangian is f + lambda.g, while Uno uses f - lambda.g
+      this->multipliers_with_flipped_sign = -multipliers;
+
+      // compute the Hessian-vector product
+      (this->asl->p.Hvcomp)(this->asl, result.data(), const_cast<double*>(x.data()), objective_number, &objective_multiplier,
+            const_cast<double*>(this->multipliers_with_flipped_sign.data()));
+      if (error_flag) {
+         throw HessianEvaluationError();
+      }
+   }
+
    double AMPLModel::variable_lower_bound(size_t variable_index) const {
       return this->variable_lower_bounds[variable_index];
    }

bindings/AMPL/AMPLModel.hpp

@@ -38,6 +38,8 @@ namespace uno {
       void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
             SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) 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/subproblems/interior_point_methods/PrimalDualInteriorPointProblem.cpp

@@ -68,6 +68,26 @@ namespace uno {
       }
    }
 
+   void PrimalDualInteriorPointProblem::compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers,
+         Vector<double>& result) const {
+      // original Lagrangian Hessian
+      this->problem.compute_hessian_vector_product(x, multipliers, result);
+
+      // barrier terms
+      for (size_t variable_index: Range(this->problem.number_variables)) {
+         double diagonal_barrier_term = 0.;
+         if (is_finite(this->problem.variable_lower_bound(variable_index))) { // lower bounded
+            const double distance_to_bound = x[variable_index] - this->problem.variable_lower_bound(variable_index);
+            diagonal_barrier_term += this->current_multipliers.lower_bounds[variable_index] / distance_to_bound;
+         }
+         if (is_finite(this->problem.variable_upper_bound(variable_index))) { // upper bounded
+            const double distance_to_bound = x[variable_index] - this->problem.variable_upper_bound(variable_index);
+            diagonal_barrier_term += this->current_multipliers.upper_bounds[variable_index] / distance_to_bound;
+         }
+         result[variable_index] += diagonal_barrier_term * x[variable_index];
+      }
+   }
+
    double PrimalDualInteriorPointProblem::variable_lower_bound(size_t /*variable_index*/) const {
       return -INF<double>;
    }

uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointProblem.hpp

@@ -19,6 +19,7 @@ namespace uno {
       void evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, const Vector<double>& multipliers,
             SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers, Vector<double>& result) 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/model/BoundRelaxedModel.hpp

@@ -32,6 +32,10 @@ namespace uno {
             SymmetricMatrix<size_t, double>& hessian) const override {
          this->model->evaluate_lagrangian_hessian(x, objective_multiplier, multipliers, hessian);
       }
+      void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) const override {
+         this->model->compute_hessian_vector_product(x, objective_multiplier, multipliers, result);
+      }
 
       // only these two functions are redefined
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override;

uno/model/FixedBoundsConstraintsModel.cpp

@@ -74,6 +74,11 @@ namespace uno {
       this->model->evaluate_lagrangian_hessian(x, objective_multiplier, multipliers, hessian);
    }
 
+   void FixedBoundsConstraintsModel::compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier,
+         const Vector<double>& multipliers, Vector<double>& result) const {
+      this->model->compute_hessian_vector_product(x, objective_multiplier, multipliers, result);
+   }
+
    double FixedBoundsConstraintsModel::variable_lower_bound(size_t variable_index) const {
       if (this->model->variable_lower_bound(variable_index) == this->model->variable_upper_bound(variable_index)) {
       // remove bounds of fixed variables

uno/model/FixedBoundsConstraintsModel.hpp

@@ -29,6 +29,8 @@ namespace uno {
       void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
             SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) 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/model/HomogeneousEqualityConstrainedModel.cpp

@@ -101,6 +101,11 @@ namespace uno {
       }
    }
 
+   void HomogeneousEqualityConstrainedModel::compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier,
+         const Vector<double>& multipliers, Vector<double>& result) const {
+      this->model->compute_hessian_vector_product(x, objective_multiplier, multipliers, result);
+   }
+
    double HomogeneousEqualityConstrainedModel::variable_lower_bound(size_t variable_index) const {
       if (variable_index < this->model->number_variables) { // original variable
          return this->model->variable_lower_bound(variable_index);

uno/model/HomogeneousEqualityConstrainedModel.hpp

@@ -26,6 +26,8 @@ namespace uno {
       void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
             SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) 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/model/Model.hpp

@@ -51,6 +51,8 @@ namespace uno {
       virtual void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const = 0;
       virtual void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
             SymmetricMatrix<size_t, double>& hessian) const = 0;
+      virtual void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) const = 0;
 
       // purely virtual functions
       [[nodiscard]] virtual double variable_lower_bound(size_t variable_index) const = 0;

uno/model/ScaledModel.cpp

@@ -11,7 +11,8 @@ namespace uno {
          Model(original_model->name + " -> scaled", original_model->number_variables, original_model->number_constraints,
                original_model->objective_sign),
          model(std::move(original_model)),
-         scaling(this->model->number_constraints, options.get_double("function_scaling_threshold")) {
+         scaling(this->model->number_constraints, options.get_double("function_scaling_threshold")),
+         scaled_multipliers(this->number_constraints) {
       if (options.get_bool("scale_functions")) {
          // evaluate the gradients at the current point
          initial_iterate.evaluate_objective_gradient(*this->model);
@@ -66,12 +67,20 @@ namespace uno {
          SymmetricMatrix<size_t, double>& hessian) const {
       // scale the objective and constraint multipliers
       const double scaled_objective_multiplier = objective_multiplier*this->scaling.get_objective_scaling();
-      // TODO preallocate this vector
-      static Vector<double> scaled_multipliers(this->number_constraints);
       for (size_t constraint_index: Range(this->number_constraints)) {
-         scaled_multipliers[constraint_index] = this->scaling.get_constraint_scaling(constraint_index) * multipliers[constraint_index];
+         this->scaled_multipliers[constraint_index] = this->scaling.get_constraint_scaling(constraint_index) * multipliers[constraint_index];
       }
-      this->model->evaluate_lagrangian_hessian(x, scaled_objective_multiplier, scaled_multipliers, hessian);
+      this->model->evaluate_lagrangian_hessian(x, scaled_objective_multiplier, this->scaled_multipliers, hessian);
+   }
+
+   void ScaledModel::compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+         Vector<double>& result) const {
+      // scale the objective and constraint multipliers
+      const double scaled_objective_multiplier = objective_multiplier*this->scaling.get_objective_scaling();
+      for (size_t constraint_index: Range(this->number_constraints)) {
+         this->scaled_multipliers[constraint_index] = this->scaling.get_constraint_scaling(constraint_index) * multipliers[constraint_index];
+      }
+      this->model->compute_hessian_vector_product(x, scaled_objective_multiplier, this->scaled_multipliers, result);
    }
 
    double ScaledModel::variable_lower_bound(size_t variable_index) const {

uno/model/ScaledModel.hpp

@@ -6,6 +6,7 @@
 
 #include <memory>
 #include "Model.hpp"
+#include "linear_algebra/Vector.hpp"
 #include "preprocessing/Scaling.hpp"
 
 namespace uno {
@@ -23,6 +24,8 @@ namespace uno {
       void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
             SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) const override;
 
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override;
       [[nodiscard]] double variable_upper_bound(size_t variable_index) const override;
@@ -53,6 +56,7 @@ namespace uno {
    private:
       const std::unique_ptr<Model> model{};
       Scaling scaling;
+      mutable Vector<double> scaled_multipliers{};
    };
 } // namespace
 

uno/optimization/EvaluationErrors.hpp

@@ -9,15 +9,21 @@ namespace uno {
       [[nodiscard]] const char* what() const noexcept override = 0;
    };
 
+   struct FunctionEvaluationError : EvaluationError {
+      [[nodiscard]] const char* what() const noexcept override {
+         return "A numerical error was encountered while evaluating a function\n";
+      }
+   };
+
    struct GradientEvaluationError : EvaluationError {
       [[nodiscard]] const char* what() const noexcept override {
          return "A numerical error was encountered while evaluating a gradient\n";
       }
    };
 
-   struct FunctionEvaluationError : EvaluationError {
+   struct HessianEvaluationError : EvaluationError {
       [[nodiscard]] const char* what() const noexcept override {
-         return "A numerical error was encountered while evaluating a function\n";
+         return "A numerical error was encountered while evaluating the Lagrangian Hessian\n";
       }
    };
 } // namespace

uno/reformulation/OptimalityProblem.cpp

@@ -32,6 +32,10 @@ namespace uno {
       this->model.evaluate_lagrangian_hessian(x, this->get_objective_multiplier(), multipliers, hessian);
    }
 
+   void OptimalityProblem::compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers, Vector<double>& result) const {
+      this->model.compute_hessian_vector_product(x, this->get_objective_multiplier(), multipliers, result);
+   }
+
    // Lagrangian gradient split in two parts: objective contribution and constraints' contribution
    void OptimalityProblem::evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient, Iterate& iterate,
          const Multipliers& multipliers) const {

uno/reformulation/OptimalityProblem.hpp

@@ -16,6 +16,7 @@ namespace uno {
       void evaluate_constraints(Iterate& iterate, std::vector<double>& constraints) const override;
       void evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, const Vector<double>& multipliers, SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers, Vector<double>& result) const override;
 
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override { return this->model.variable_lower_bound(variable_index); }
       [[nodiscard]] double variable_upper_bound(size_t variable_index) const override { return this->model.variable_upper_bound(variable_index); }

uno/reformulation/OptimizationProblem.hpp

@@ -42,6 +42,7 @@ namespace uno {
       virtual void evaluate_constraints(Iterate& iterate, std::vector<double>& constraints) const = 0;
       virtual void evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix<double>& constraint_jacobian) const = 0;
       virtual void evaluate_lagrangian_hessian(const Vector<double>& x, const Vector<double>& multipliers, SymmetricMatrix<size_t, double>& hessian) const = 0;
+      virtual void compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers, Vector<double>& result) const = 0;
 
       [[nodiscard]] size_t get_number_original_variables() const;
       [[nodiscard]] virtual double variable_lower_bound(size_t variable_index) const = 0;

uno/reformulation/l1RelaxedProblem.cpp

@@ -113,6 +113,19 @@ namespace uno {
       }
    }
 
+   void l1RelaxedProblem::compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers, Vector<double>& result) const {
+      this->model.compute_hessian_vector_product(x, this->objective_multiplier, multipliers, result);
+
+      // proximal contribution
+      if (this->proximal_center != nullptr && this->proximal_coefficient != 0.) {
+         for (size_t variable_index: Range(this->model.number_variables)) {
+            const double scaling = std::min(1., 1./std::abs(this->proximal_center[variable_index]));
+            const double proximal_term = this->proximal_coefficient * scaling * scaling;
+            result[variable_index] += proximal_term * x[variable_index];
+         }
+      }
+   }
+
    // Lagrangian gradient split in two parts: objective contribution and constraints' contribution
    void l1RelaxedProblem::evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient, Iterate& iterate,
          const Multipliers& multipliers) const {

uno/reformulation/l1RelaxedProblem.hpp

@@ -19,6 +19,7 @@ namespace uno {
       void evaluate_constraints(Iterate& iterate, std::vector<double>& constraints) const override;
       void evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, const Vector<double>& multipliers, SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers, Vector<double>& result) 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/solvers/BQPD/BQPDSolver.cpp

@@ -20,7 +20,7 @@
 #define hessian_vector_product FC_GLOBAL(gdotx, GDOTX)
 
 extern "C" {
-   void hessian_vector_product(int *n, const double x[], const double ws[], const int lws[], double v[]);
+   void hessian_vector_product([[maybe_unused]] int *n, const double x[], const double ws[], const int lws[], double v[]);
 
    // fortran common block used in bqpd/bqpd.f
    extern struct {
@@ -107,7 +107,7 @@ namespace uno {
       WSC.ll = static_cast<int>(this->size_hessian_sparsity);
       WSC.mxws = static_cast<int>(this->size_hessian_workspace);
       WSC.mxlws = static_cast<int>(this->size_hessian_sparsity_workspace);
-      ALPHAC.alpha = 0; // inertia control
+      ALPHAC.alpha = 0.; // inertia control
 
       if (this->print_subproblem) {
          DEBUG << "objective gradient: " << linear_objective;