[WIP] L-BFGS approximation

.github/julia/build_tarballs_release.jl

@@ -167,4 +167,4 @@ build_tarballs(
     julia_compat = "1.6",
     preferred_gcc_version = v"10.2.0",
     clang_use_lld=false,
-)
+)

.github/julia/build_tarballs_yggdrasil.jl

@@ -112,4 +112,4 @@ build_tarballs(
     julia_compat = "1.9",
     preferred_gcc_version = v"10.2.0",
     clang_use_lld=false,
-)
+)

CMakeLists.txt

@@ -44,7 +44,9 @@ find_library(HSL hsl)
 find_library(MA57 ma57)
 find_library(MA27 ma27)
 find_library(BQPD bqpd)
-if(HSL OR MA57 OR MA27 OR BQPD)
+option(WITH_LAPACK "Enable LAPACK" OFF)
+find_package(LAPACK)
+if(HSL OR MA57 OR MA27 OR BQPD OR (WITH_LAPACK AND LAPACK_FOUND))
    message(STATUS "Fortran compiler required")
    enable_language(Fortran)
    include(FortranCInterface)
@@ -89,6 +91,7 @@ file(GLOB TESTS_UNO_SOURCE_FILES
    unotest/unit_tests/ConcatenationTests.cpp
    unotest/unit_tests/COOSparseStorageTests.cpp
    unotest/unit_tests/CSCSparseStorageTests.cpp
+   unotest/unit_tests/DenseMatrixTests.cpp
    unotest/unit_tests/RangeTests.cpp
    unotest/unit_tests/ScalarMultipleTests.cpp
    unotest/unit_tests/SparseVectorTests.cpp
@@ -101,6 +104,7 @@ file(GLOB TESTS_UNO_SOURCE_FILES
 # external dependencies #
 #########################
 set(LIBRARIES "")
+set(TEST_LIBRARIES "")
 
 # function that links an existing library to Uno
 function(link_to_uno library_name library_path)
@@ -144,6 +148,20 @@ if(HSL OR MA27)
    list(APPEND LIBRARIES ${BLAS_LIBRARIES})
 endif()
 
+# LAPACK
+if(WITH_LAPACK)
+	find_package(LAPACK)
+	if(NOT LAPACK_FOUND)
+		message(WARNING "Optional library LAPACK was not found.")
+	else()
+		link_to_uno(LAPACK ${LAPACK_LIBRARIES})
+		set(TEST_LIBRARIES ${TEST_LIBRARIES} ${LAPACK_LIBRARIES})
+		# enable the use of quasi-Newton methods
+		list(APPEND UNO_SOURCE_FILES uno/ingredients/hessian_models/quasi_newton/LBFGSHessian.cpp)
+		list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/LAPACKTests.cpp)
+	endif()
+endif()
+
 # METIS
 find_library(METIS metis)
 if(METIS)

uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp

@@ -30,8 +30,8 @@ namespace uno {
          optimality_problem(model),
          // relax the linear constraints in the l1 relaxed problem only if we are using a trust-region constraint
          feasibility_problem(model, 0., this->constraint_violation_coefficient, use_trust_region),
-         optimality_hessian_model(HessianModelFactory::create(model, options)),
-         feasibility_hessian_model(HessianModelFactory::create(model, options)),
+         optimality_hessian_model(HessianModelFactory::create(model, 1., options)),
+         feasibility_hessian_model(HessianModelFactory::create(model, 0., options)),
          optimality_inertia_correction_strategy(InertiaCorrectionStrategyFactory::create(options)),
          feasibility_inertia_correction_strategy(InertiaCorrectionStrategyFactory::create(options)),
          optimality_inequality_handling_method(InequalityHandlingMethodFactory::create(options)),
@@ -59,6 +59,8 @@ namespace uno {
       this->optimality_inequality_handling_method->initialize_statistics(statistics);
       this->feasibility_inertia_correction_strategy->initialize_statistics(statistics);
       this->feasibility_inequality_handling_method->initialize_statistics(statistics);
+      this->optimality_hessian_model->initialize_statistics(statistics);
+      this->feasibility_hessian_model->initialize_statistics(statistics);
       statistics.add_column("Phase", Statistics::int_width - 1, 3, Statistics::column_order.at("Phase"));
       statistics.set("Phase", "OPT");
 
@@ -201,10 +203,16 @@ namespace uno {
       if (this->current_phase == Phase::OPTIMALITY) {
          accept_iterate = this->optimality_inequality_handling_method->is_iterate_acceptable(statistics,
             *this->optimality_globalization_strategy, current_iterate, trial_iterate, direction, step_length, user_callbacks);
+         if (accept_iterate) {
+            this->optimality_hessian_model->notify_accepted_iterate(statistics, current_iterate, trial_iterate);
+         }
       }
       else {
          accept_iterate = this->feasibility_inequality_handling_method->is_iterate_acceptable(statistics,
             this->feasibility_globalization_strategy, current_iterate, trial_iterate, direction, step_length, user_callbacks);
+         if (accept_iterate) {
+            this->feasibility_hessian_model->notify_accepted_iterate(statistics, current_iterate, trial_iterate);
+         }
       }
       trial_iterate.status = this->check_termination(model, trial_iterate);
 

uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp

@@ -17,7 +17,7 @@ namespace uno {
          ConstraintRelaxationStrategy(options),
          problem(model),
          inequality_handling_method(InequalityHandlingMethodFactory::create(options)),
-         hessian_model(HessianModelFactory::create(model, options)),
+         hessian_model(HessianModelFactory::create(model, 1., options)),
          inertia_correction_strategy(InertiaCorrectionStrategyFactory::create(options)),
          globalization_strategy(options) {
    }
@@ -32,6 +32,7 @@ namespace uno {
       // statistics
       this->inertia_correction_strategy->initialize_statistics(statistics);
       this->inequality_handling_method->initialize_statistics(statistics);
+      this->hessian_model->initialize_statistics(statistics);
 
       // initial iterate
       this->inequality_handling_method->generate_initial_iterate(initial_iterate);
@@ -69,6 +70,9 @@ namespace uno {
          UserCallbacks& user_callbacks) {
       const bool accept_iterate = this->inequality_handling_method->is_iterate_acceptable(statistics, this->globalization_strategy,
          current_iterate, trial_iterate, direction, step_length, user_callbacks);
+      if (accept_iterate) {
+         this->hessian_model->notify_accepted_iterate(statistics, current_iterate, trial_iterate);
+      }
       trial_iterate.status = this->check_termination(model, trial_iterate);
       warmstart_information.no_changes();
       return accept_iterate;

uno/ingredients/globalization_strategies/switching_methods/filter_methods/WaechterFilterMethod.cpp

@@ -90,4 +90,4 @@ namespace uno {
    std::string WaechterFilterMethod::get_name() const {
       return "Waechter-filter";
    }
-} // namespace
+} // namespace

uno/ingredients/hessian_models/ExactHessian.cpp

@@ -33,6 +33,12 @@ namespace uno {
       return false;
    }
 
+   void ExactHessian::initialize_statistics(Statistics& /*statistics*/) const {
+   }
+
+   void ExactHessian::notify_accepted_iterate(Statistics& /*statistics*/, Iterate& /*current_iterate*/, Iterate& /*trial_iterate*/) {
+   }
+
    void ExactHessian::evaluate_hessian(Statistics& /*statistics*/, const Vector<double>& primal_variables,
          double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values) {
       this->model.evaluate_lagrangian_hessian(primal_variables, objective_multiplier, constraint_multipliers, hessian_values);

uno/ingredients/hessian_models/ExactHessian.hpp

@@ -13,6 +13,7 @@ namespace uno {
    class ExactHessian : public HessianModel {
    public:
       explicit ExactHessian(const Model& model);
+      ~ExactHessian() override = default;
 
       [[nodiscard]] bool has_hessian_operator() const override;
       [[nodiscard]] bool has_hessian_matrix() const override;
@@ -21,6 +22,8 @@ namespace uno {
       void compute_sparsity(uno_int* row_indices, uno_int* column_indices, uno_int solver_indexing) const override;
       [[nodiscard]] bool is_positive_definite() const override;
 
+      void initialize_statistics(Statistics& statistics) const override;
+      void notify_accepted_iterate(Statistics& statistics, Iterate& current_iterate, Iterate& trial_iterate) override;
       void evaluate_hessian(Statistics& statistics, const Vector<double>& primal_variables, double objective_multiplier,
          const Vector<double>& constraint_multipliers, double* hessian_values) override;
       void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier,

uno/ingredients/hessian_models/HessianModel.hpp

@@ -11,6 +11,9 @@
 
 namespace uno {
    // forward declarations
+   class Iterate;
+   class Model;
+   class Options;
    class Statistics;
    template <typename ElementType>
    class Vector;
@@ -31,6 +34,8 @@ namespace uno {
       virtual void compute_sparsity(uno_int* row_indices, uno_int* column_indices, uno_int solver_indexing) const = 0;
       [[nodiscard]] virtual bool is_positive_definite() const = 0;
 
+      virtual void initialize_statistics(Statistics& statistics) const = 0;
+      virtual void notify_accepted_iterate(Statistics& statistics, Iterate& current_iterate, Iterate& trial_iterate) = 0;
       virtual void evaluate_hessian(Statistics& statistics, const Vector<double>& primal_variables,
          double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values) = 0;
       virtual void compute_hessian_vector_product(const double* x, const double* vector,

uno/ingredients/hessian_models/HessianModelFactory.cpp

@@ -8,12 +8,16 @@
 #include "ExactHessian.hpp"
 #include "IdentityHessian.hpp"
 #include "ZeroHessian.hpp"
+#ifdef HAS_LAPACK
+#include "quasi_newton/LBFGSHessian.hpp"
+#endif
 #include "model/Model.hpp"
 #include "options/Options.hpp"
 #include "tools/Logger.hpp"
 
 namespace uno {
-   std::unique_ptr<HessianModel> HessianModelFactory::create(const Model& model, const Options& options) {
+   std::unique_ptr<HessianModel> HessianModelFactory::create(const Model& model, [[maybe_unused]] double objective_multiplier,
+         const Options& options) {
       const std::string& hessian_model = options.get_string("hessian_model");
       if (hessian_model == "exact") {
          // if no Hessian (matrix or operator) is available, pick a zero Hessian
@@ -28,6 +32,11 @@ namespace uno {
       else if (hessian_model == "identity") {
          return std::make_unique<IdentityHessian>(model.number_variables);
       }
+#ifdef HAS_LAPACK
+      else if (hessian_model == "L-BFGS") {
+         return std::make_unique<LBFGSHessian>(model, objective_multiplier, options);
+      }
+#endif
       else if (hessian_model == "zero") {
          return std::make_unique<ZeroHessian>(model.number_variables);
       }

uno/ingredients/hessian_models/HessianModelFactory.hpp

@@ -14,7 +14,7 @@ namespace uno {
 
    class HessianModelFactory {
    public:
-      static std::unique_ptr<HessianModel> create(const Model& model, const Options& options);
+      static std::unique_ptr<HessianModel> create(const Model& model, double objective_multiplier, const Options& options);
    };
 } // namespace
 

uno/ingredients/hessian_models/IdentityHessian.cpp

@@ -37,6 +37,12 @@ namespace uno {
       return true;
    }
 
+   void IdentityHessian::initialize_statistics(Statistics& /*statistics*/) const {
+   }
+
+   void IdentityHessian::notify_accepted_iterate(Statistics& /*statistics*/, Iterate& /*current_iterate*/, Iterate& /*trial_iterate*/) {
+   }
+
    void IdentityHessian::evaluate_hessian(Statistics& /*statistics*/, const Vector<double>& /*primal_variables*/,
          double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* hessian_values) {
       DEBUG << "Setting identity Hessian\n";

uno/ingredients/hessian_models/IdentityHessian.hpp

@@ -18,6 +18,8 @@ namespace uno {
       void compute_sparsity(uno_int* row_indices, uno_int* column_indices, uno_int solver_indexing) const override;
       [[nodiscard]] bool is_positive_definite() const override;
 
+      void initialize_statistics(Statistics& statistics) const override;
+      void notify_accepted_iterate(Statistics& statistics, Iterate& current_iterate, Iterate& trial_iterate) override;
       void evaluate_hessian(Statistics& statistics, const Vector<double>& primal_variables,
          double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values) override;
       void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier,

uno/ingredients/hessian_models/ZeroHessian.cpp

@@ -32,9 +32,14 @@ namespace uno {
       return false;
    }
 
+   void ZeroHessian::initialize_statistics(Statistics& /*statistics*/) const {
+   }
+
+   void ZeroHessian::notify_accepted_iterate(Statistics& /*statistics*/, Iterate& /*current_iterate*/, Iterate& /*trial_iterate*/) {
+   }
+
    void ZeroHessian::evaluate_hessian(Statistics& /*statistics*/, const Vector<double>& /*primal_variables*/,
          double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* /*hessian_values*/) {
-      // do nothing
    }
 
    void ZeroHessian::compute_hessian_vector_product(const double* /*x*/, const double* /*vector*/,

uno/ingredients/hessian_models/ZeroHessian.hpp

@@ -10,6 +10,7 @@ namespace uno {
    class ZeroHessian : public HessianModel {
    public:
       explicit ZeroHessian(size_t number_variables);
+      ~ZeroHessian() override = default;
 
       [[nodiscard]] bool has_hessian_operator() const override;
       [[nodiscard]] bool has_hessian_matrix() const override;
@@ -18,6 +19,8 @@ namespace uno {
       void compute_sparsity(uno_int* row_indices, uno_int* column_indices, uno_int solver_indexing) const override;
       [[nodiscard]] bool is_positive_definite() const override;
 
+      void initialize_statistics(Statistics& statistics) const override;
+      void notify_accepted_iterate(Statistics& statistics, Iterate& current_iterate, Iterate& trial_iterate) override;
       void evaluate_hessian(Statistics& statistics, const Vector<double>& primal_variables, double objective_multiplier,
          const Vector<double>& constraint_multipliers, double* hessian_values) override;
       void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier,