Some progress on L-BFGS

cvanaret · cvanaret · commit b3cc74e3ccb8 · 2025-04-05T13:36:36.000+02:00
diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp
@@ -178,7 +178,7 @@ namespace uno {
       }
       ConstraintRelaxationStrategy::set_progress_statistics(statistics, trial_iterate);
       if (accept_iterate) {
-         this->inequality_handling_method->notify_accepted_iterate(trial_iterate);
+         this->inequality_handling_method->notify_accepted_iterate(current_iterate, trial_iterate);
          user_callbacks.notify_acceptable_iterate(trial_iterate.primals,
                this->current_phase == Phase::OPTIMALITY ? trial_iterate.multipliers : trial_iterate.feasibility_multipliers,
                this->current_problem().get_objective_multiplier());
diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
@@ -262,7 +262,7 @@ namespace uno {
       if (accept_iterate) {
          this->check_exact_relaxation(trial_iterate);
          // this->set_dual_residuals_statistics(statistics, trial_iterate);
-         this->inequality_handling_method->notify_accepted_iterate(trial_iterate);
+         this->inequality_handling_method->notify_accepted_iterate(current_iterate, trial_iterate);
          user_callbacks.notify_acceptable_iterate(trial_iterate.primals, trial_iterate.multipliers, this->penalty_parameter);
       }
       this->set_progress_statistics(statistics, trial_iterate);
diff --git a/uno/ingredients/hessian_models/ConvexifiedHessian.cpp b/uno/ingredients/hessian_models/ConvexifiedHessian.cpp
@@ -25,7 +25,7 @@ namespace uno {
       statistics.add_column("regulariz", Statistics::double_width - 4, options.get_int("statistics_regularization_column_order"));
    }
 
-   void ConvexifiedHessian::notify_accepted_iterate(const Iterate& /*iterate*/) {
+   void ConvexifiedHessian::notify_accepted_iterate(const Iterate& /*current_iterate*/, const Iterate& /*trial_iterate*/) {
       // do nothing
    }
 
diff --git a/uno/ingredients/hessian_models/ConvexifiedHessian.hpp b/uno/ingredients/hessian_models/ConvexifiedHessian.hpp
@@ -16,7 +16,7 @@ namespace uno {
       ConvexifiedHessian(size_t dimension, size_t number_nonzeros, const Options& options);
 
       void initialize_statistics(Statistics& statistics, const Options& options) const override;
-      void notify_accepted_iterate(const Iterate& iterate) override;
+      void notify_accepted_iterate(const Iterate& current_iterate, const Iterate& trial_iterate) override;
       void evaluate_hessian(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
             const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) override;
       void compute_hessian_vector_product(const OptimizationProblem& problem, const Vector<double>& vector,
diff --git a/uno/ingredients/hessian_models/ExactHessian.cpp b/uno/ingredients/hessian_models/ExactHessian.cpp
@@ -14,7 +14,7 @@ namespace uno {
       // do nothing
    }
 
-   void ExactHessian::notify_accepted_iterate(const Iterate& /*iterate*/) {
+   void ExactHessian::notify_accepted_iterate(const Iterate& /*current_iterate*/, const Iterate& /*trial_iterate*/) {
       // do nothing
    }
 
diff --git a/uno/ingredients/hessian_models/ExactHessian.hpp b/uno/ingredients/hessian_models/ExactHessian.hpp
@@ -10,7 +10,7 @@ namespace uno {
       ExactHessian();
 
       void initialize_statistics(Statistics& statistics, const Options& options) const override;
-      void notify_accepted_iterate(const Iterate& iterate) override;
+      void notify_accepted_iterate(const Iterate& current_iterate, const Iterate& trial_iterate) override;
       void evaluate_hessian(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
             const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) override;
       void compute_hessian_vector_product(const OptimizationProblem& problem, const Vector<double>& vector,
diff --git a/uno/ingredients/hessian_models/HessianModel.hpp b/uno/ingredients/hessian_models/HessianModel.hpp
@@ -25,7 +25,7 @@ namespace uno {
       size_t evaluation_count{0};
 
       virtual void initialize_statistics(Statistics& statistics, const Options& options) const = 0;
-      virtual void notify_accepted_iterate(const Iterate& iterate) = 0;
+      virtual void notify_accepted_iterate(const Iterate& current_iterate, const Iterate& trial_iterate) = 0;
       virtual void evaluate_hessian(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
             const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) = 0;
       virtual void compute_hessian_vector_product(const OptimizationProblem& problem, const Vector<double>& vector,
diff --git a/uno/ingredients/hessian_models/LBFGSHessian.cpp b/uno/ingredients/hessian_models/LBFGSHessian.cpp
@@ -3,6 +3,8 @@
 
 #include <stdexcept>
 #include "LBFGSHessian.hpp"
+#include "optimization/Iterate.hpp"
+#include "symbolic/Range.hpp"
 
 #ifdef WITH_LAPACK
 #include "fortran_interface.h"
@@ -16,18 +18,50 @@ namespace uno {
    LBFGSHessian::LBFGSHessian(size_t dimension, size_t memory_size):
          HessianModel(),
          dimension(dimension),
-         memory_size(memory_size),
+         memory_capacity(memory_size),
          S_matrix(dimension, memory_size),
          Y_matrix(dimension, memory_size),
+         L_matrix(memory_size, memory_size),
+         D_matrix(memory_size),
          M_matrix(memory_size, memory_size) {
    }
 
    void LBFGSHessian::initialize_statistics(Statistics& /*statistics*/, const Options& /*options*/) const {
       // do nothing
    }
 
-   void LBFGSHessian::notify_accepted_iterate(const Iterate& /*iterate*/) {
+   void LBFGSHessian::notify_accepted_iterate(const Iterate& current_iterate, const Iterate& trial_iterate) {
+      std::cout << "Adding vector to L-BFGS memory at slot " << this->current_available_slot << '\n';
+      // this->current_available_slot lives in [0, this->memory_capacity)
 
+      // TODO figure out if we're extending or replacing in memory
+
+      // fill the S matrix
+      for (size_t variable_index: Range(this->dimension)) {
+         this->S_matrix.entry(variable_index, this->current_available_slot) = trial_iterate.primals[variable_index] -
+            current_iterate.primals[variable_index];
+      }
+      std::cout << "S:\n" << this->S_matrix;
+
+      // fill the Y matrix
+      // TODO
+
+      // fill the D matrix (diagonal)
+      this->D_matrix[this->current_available_slot] = 1.; // TODO dot(s_new, y_new)
+
+      // fill the L matrix (lower triangular with a zero diagonal)
+      for (size_t column_index: Range(this->current_memory_size)) {
+         for (size_t row_index: Range(column_index+1, this->current_memory_size)) {
+            this->L_matrix.entry(row_index, column_index) = 1.; // TODO dot(s_i, y_j)
+         }
+      }
+      std::cout << "L:\n" << this->L_matrix;
+
+      // if we exceed the size of the memory, we start over and replace the older point in memory
+      this->current_available_slot = (this->current_available_slot + 1) % this->memory_capacity;
+      this->current_memory_size = std::min(current_memory_size + 1, this->memory_capacity);
+      std::cout << "There are now " << this->current_memory_size << " iterates in memory (capacity " <<
+         this->memory_capacity << ")\n";
    }
 
    void LBFGSHessian::evaluate_hessian(Statistics& /*statistics*/, const OptimizationProblem& /*problem*/,
@@ -36,8 +70,11 @@ namespace uno {
       throw std::runtime_error("LBFGSHessian::evaluate_hessian not implemented");
    }
 
-   void LBFGSHessian::compute_hessian_vector_product(const OptimizationProblem& /*problem*/, const Vector<double>& /*vector*/,
-         const Vector<double>& /*constraint_multipliers*/, Vector<double>& /*result*/) {
-      throw std::runtime_error("LBFGSHessian::compute_hessian_vector_product not implemented");
+   void LBFGSHessian::compute_hessian_vector_product(const OptimizationProblem& /*problem*/, const Vector<double>& vector,
+         const Vector<double>& /*constraint_multipliers*/, Vector<double>& result) {
+      //throw std::runtime_error("LBFGSHessian::compute_hessian_vector_product not implemented");
+      for (size_t variable_index: Range(vector.size())) {
+         result[variable_index] = vector[variable_index];
+      }
    }
 } // namespace
diff --git a/uno/ingredients/hessian_models/LBFGSHessian.hpp b/uno/ingredients/hessian_models/LBFGSHessian.hpp
@@ -14,17 +14,21 @@ namespace uno {
       ~LBFGSHessian() override = default;
 
       void initialize_statistics(Statistics& statistics, const Options& options) const override;
-      void notify_accepted_iterate(const Iterate& iterate) override;
+      void notify_accepted_iterate(const Iterate& current_iterate, const Iterate& trial_iterate) override;
       void evaluate_hessian(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
             const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) override;
       void compute_hessian_vector_product(const OptimizationProblem& problem, const Vector<double>& vector,
          const Vector<double>& constraint_multipliers, Vector<double>& result) override;
 
    protected:
       const size_t dimension;
-      const size_t memory_size;
+      const size_t memory_capacity;
+      size_t current_memory_size{0};
+      size_t current_available_slot{0};
       DenseMatrix<double> S_matrix;
       DenseMatrix<double> Y_matrix;
+      DenseMatrix<double> L_matrix;
+      std::vector<double> D_matrix; // D is diagonal
       DenseMatrix<double> M_matrix;
    };
 } // namespace
diff --git a/uno/ingredients/hessian_models/ZeroHessian.cpp b/uno/ingredients/hessian_models/ZeroHessian.cpp
@@ -8,7 +8,7 @@
 namespace uno {
    void ZeroHessian::initialize_statistics(Statistics& /*statistics*/, const Options& /*options*/) const { }
 
-   void ZeroHessian::notify_accepted_iterate(const Iterate& /*iterate*/) {
+   void ZeroHessian::notify_accepted_iterate(const Iterate& /*current_iterate*/, const Iterate& /*trial_iterate*/) {
       // do nothing
    }
 
diff --git a/uno/ingredients/hessian_models/ZeroHessian.hpp b/uno/ingredients/hessian_models/ZeroHessian.hpp
@@ -10,7 +10,7 @@ namespace uno {
       ZeroHessian() = default;
 
       void initialize_statistics(Statistics& statistics, const Options& options) const override;
-      void notify_accepted_iterate(const Iterate& iterate) override;
+      void notify_accepted_iterate(const Iterate& current_iterate, const Iterate& trial_iterate) override;
       void evaluate_hessian(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
             const Vector<double>& constraint_multipliers, SymmetricMatrix<size_t, double>& hessian) override;
       void compute_hessian_vector_product(const OptimizationProblem& problem, const Vector<double>& vector,
diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp
@@ -16,8 +16,8 @@ namespace uno {
       this->trust_region_radius = new_trust_region_radius;
    }
 
-   void InequalityHandlingMethod::notify_accepted_iterate(const Iterate& iterate) const {
-      this->hessian_model->notify_accepted_iterate(iterate);
+   void InequalityHandlingMethod::notify_accepted_iterate(const Iterate& current_iterate, const Iterate& trial_iterate) const {
+      this->hessian_model->notify_accepted_iterate(current_iterate, trial_iterate);
    }
 
    size_t InequalityHandlingMethod::get_hessian_evaluation_count() const {
diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp
@@ -48,7 +48,7 @@ namespace uno {
       [[nodiscard]] virtual double compute_predicted_auxiliary_reduction_model(const Model& model, const Iterate& current_iterate,
             const Vector<double>& primal_direction, double step_length) const = 0;
 
-      void notify_accepted_iterate(const Iterate& iterate) const;
+      void notify_accepted_iterate(const Iterate &current_iterate, const Iterate& trial_iterate) const;
       virtual void postprocess_iterate(const OptimizationProblem& problem, Iterate& iterate) = 0;
 
       [[nodiscard]] size_t get_hessian_evaluation_count() const;
diff --git a/uno/linear_algebra/DenseMatrix.hpp b/uno/linear_algebra/DenseMatrix.hpp
@@ -5,7 +5,9 @@
 #define UNO_DENSEMATRIX_H
 
 #include <cstddef>
+#include <iostream>
 #include <vector>
+#include "symbolic/Range.hpp"
 
 namespace uno {
    template <typename ElementType>
@@ -14,9 +16,12 @@ namespace uno {
       DenseMatrix(size_t number_rows, size_t number_columns);
       ~DenseMatrix() = default;
 
-      ElementType& get(size_t row_index, size_t column_index);
+      ElementType& entry(size_t row_index, size_t column_index);
+      const ElementType& entry(size_t row_index, size_t column_index) const;
       ElementType* data() const;
 
+      void print(std::ostream& stream) const;
+
    protected:
       const size_t number_rows, number_columns;
       std::vector<ElementType> matrix; // column-major ordering
@@ -29,14 +34,37 @@ namespace uno {
    }
 
    template <typename ElementType>
-   ElementType& DenseMatrix<ElementType>::get(size_t row_index, size_t column_index) {
+   ElementType& DenseMatrix<ElementType>::entry(size_t row_index, size_t column_index) {
+      return this->matrix[column_index * this->number_rows + row_index];
+   }
+
+   template <typename ElementType>
+   const ElementType& DenseMatrix<ElementType>::entry(size_t row_index, size_t column_index) const {
       return this->matrix[column_index * this->number_rows + row_index];
    }
 
    template <typename ElementType>
    ElementType* DenseMatrix<ElementType>::data() const {
       return this->matrix.data();
    }
+
+   template <typename ElementType>
+   void DenseMatrix<ElementType>::print(std::ostream& stream) const {
+      stream << "Dense matrix (" << this->number_rows << "x" << this->number_columns << ")\n";
+      for (size_t column_index: Range(this->number_columns)) {
+         stream << "Column " << column_index << ":";
+         for (size_t row_index: Range(this->number_rows)) {
+            stream << ' ' << this->entry(row_index, column_index);
+         }
+         stream << '\n';
+      }
+   }
+
+   template <typename ElementType>
+   std::ostream& operator<<(std::ostream& stream, const DenseMatrix<ElementType>& matrix) {
+      matrix.print(stream);
+      return stream;
+   }
 } // namespace
 
 #endif // UNO_DENSEMATRIX_H
diff --git a/unotest/unit_tests/DenseMatrixTests.cpp b/unotest/unit_tests/DenseMatrixTests.cpp
@@ -8,5 +8,5 @@ using namespace uno;
 
 TEST(DenseMatrix, Access) {
    DenseMatrix<double> M(5, 2);
-   ASSERT_EQ(M.get(1, 1), 0.);
+   ASSERT_EQ(M.entry(1, 1), 0.);
 }

Original file line number	Diff line number	Diff line change
`@@ -178,7 +178,7 @@ namespace uno {`
`178`	`178`	`}`
`179`	`179`	`ConstraintRelaxationStrategy::set_progress_statistics(statistics, trial_iterate);`
`180`	`180`	`if (accept_iterate) {`
`181`		`- this->inequality_handling_method->notify_accepted_iterate(trial_iterate);`
	`181`	`+ this->inequality_handling_method->notify_accepted_iterate(current_iterate, trial_iterate);`
`182`	`182`	`user_callbacks.notify_acceptable_iterate(trial_iterate.primals,`
`183`	`183`	`this->current_phase == Phase::OPTIMALITY ? trial_iterate.multipliers : trial_iterate.feasibility_multipliers,`
`184`	`184`	`this->current_problem().get_objective_multiplier());`
Original file line number	Diff line number	Diff line change
`@@ -262,7 +262,7 @@ namespace uno {`
`262`	`262`	`if (accept_iterate) {`
`263`	`263`	`this->check_exact_relaxation(trial_iterate);`
`264`	`264`	`// this->set_dual_residuals_statistics(statistics, trial_iterate);`
`265`		`- this->inequality_handling_method->notify_accepted_iterate(trial_iterate);`
	`265`	`+ this->inequality_handling_method->notify_accepted_iterate(current_iterate, trial_iterate);`
`266`	`266`	`user_callbacks.notify_acceptable_iterate(trial_iterate.primals, trial_iterate.multipliers, this->penalty_parameter);`
`267`	`267`	`}`
`268`	`268`	`this->set_progress_statistics(statistics, trial_iterate);`
Original file line number	Diff line number	Diff line change
`@@ -25,7 +25,7 @@ namespace uno {`
`25`	`25`	`statistics.add_column("regulariz", Statistics::double_width - 4, options.get_int("statistics_regularization_column_order"));`
`26`	`26`	`}`
`27`	`27`
`28`		`- void ConvexifiedHessian::notify_accepted_iterate(const Iterate& /iterate/) {`
	`28`	`+ void ConvexifiedHessian::notify_accepted_iterate(const Iterate& /current_iterate/, const Iterate& /trial_iterate/) {`
`29`	`29`	`// do nothing`
`30`	`30`	`}`
`31`	`31`
Original file line number	Diff line number	Diff line change
`@@ -14,7 +14,7 @@ namespace uno {`
`14`	`14`	`// do nothing`
`15`	`15`	`}`
`16`	`16`
`17`		`- void ExactHessian::notify_accepted_iterate(const Iterate& /iterate/) {`
	`17`	`+ void ExactHessian::notify_accepted_iterate(const Iterate& /current_iterate/, const Iterate& /trial_iterate/) {`
`18`	`18`	`// do nothing`
`19`	`19`	`}`
`20`	`20`
Original file line number	Diff line number	Diff line change
`@@ -8,7 +8,7 @@`
`8`	`8`	`namespace uno {`
`9`	`9`	`void ZeroHessian::initialize_statistics(Statistics& /statistics/, const Options& /options/) const { }`
`10`	`10`
`11`		`- void ZeroHessian::notify_accepted_iterate(const Iterate& /iterate/) {`
	`11`	`+ void ZeroHessian::notify_accepted_iterate(const Iterate& /current_iterate/, const Iterate& /trial_iterate/) {`
`12`	`12`	`// do nothing`
`13`	`13`	`}`
`14`	`14`