[WIP] Dogleg method

CMakeLists.txt

@@ -76,6 +76,7 @@ file(GLOB UNO_SOURCE_FILES
    uno/ingredients/inertia_correction_strategies/*.cpp
    uno/ingredients/subproblem/*.cpp
    uno/ingredients/subproblem_solvers/*.cpp
+   uno/ingredients/subproblem_solvers/dogleg/*.cpp
    uno/model/*.cpp
    uno/optimization/*.cpp
    uno/options/*.cpp

uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp

@@ -0,0 +1,32 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#ifndef UNO_TRUSTREGIONSOLVER_H
+#define UNO_TRUSTREGIONSOLVER_H
+
+#include "InequalityConstrainedSolver.hpp"
+
+namespace uno {
+   // forward declarations
+   class Direction;
+   class Statistics;
+   class Subproblem;
+   template <typename ElementType>
+   class Vector;
+   class WarmstartInformation;
+
+   class BoundConstrainedSolver: public InequalityConstrainedSolver {
+   public:
+      BoundConstrainedSolver() = default;
+      ~BoundConstrainedSolver() override = default;
+
+      void initialize_memory(const Subproblem& subproblem) override = 0;
+
+      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector<double>& initial_point,
+         Direction& direction, const WarmstartInformation& warmstart_information) override = 0;
+
+      [[nodiscard]] virtual EvaluationSpace& get_evaluation_space() = 0;
+   };
+} // namespace
+
+#endif // UNO_TRUSTREGIONSOLVER_H

uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp

@@ -0,0 +1,93 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#include "DoglegEvaluationSpace.hpp"
+#include "ingredients/subproblem/Subproblem.hpp"
+#include "ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp"
+#include "optimization/Direction.hpp"
+#include "optimization/OptimizationProblem.hpp"
+#include "optimization/WarmstartInformation.hpp"
+
+namespace uno {
+   void DoglegEvaluationSpace::initialize_memory(const Subproblem& subproblem) {
+      // Newton step
+      this->objective_gradient.resize(subproblem.number_variables);
+      this->newton_step.resize(subproblem.number_variables);
+      // Cauchy step
+      this->hessian_gradient_product.resize(subproblem.number_variables);
+      this->cauchy_step.resize(subproblem.number_variables);
+   }
+
+   void DoglegEvaluationSpace::evaluate_constraint_jacobian(const OptimizationProblem& /*problem*/, Iterate& /*iterate*/) {
+      // do nothing
+   }
+
+   void DoglegEvaluationSpace::compute_constraint_jacobian_vector_product(const Vector<double>& /*vector*/, Vector<double>& result) const {
+      result.fill(0.);
+   }
+
+   void DoglegEvaluationSpace::compute_constraint_jacobian_transposed_vector_product(const Vector<double>& /*vector*/,
+         Vector<double>& result) const {
+      result.fill(0.);
+   }
+
+   double DoglegEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/,
+         const Vector<double>& /*vector*/) const {
+      throw std::runtime_error("Not implemented yet");
+   }
+
+   void DoglegEvaluationSpace::evaluate_objective_gradient(const Subproblem& subproblem, const WarmstartInformation& warmstart_information) {
+      if (warmstart_information.objective_changed) {
+         subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->objective_gradient.data());
+      }
+   }
+
+   void DoglegEvaluationSpace::compute_newton_step(Statistics& statistics, const Subproblem& subproblem,
+         SymmetricIndefiniteLinearSolver<double>& linear_solver, Direction& direction, const WarmstartInformation& warmstart_information) {
+      if (warmstart_information.objective_changed) {
+         // g = ∇f(x_k)
+         subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->objective_gradient.data());
+         linear_solver.solve_indefinite_system(statistics, subproblem, direction, warmstart_information);
+         this->newton_step = direction.primals;
+         this->newton_step_squared_norm = dot(this->newton_step, this->newton_step);
+      }
+   }
+
+   void DoglegEvaluationSpace::compute_dogleg(const Subproblem& subproblem, Direction& /*direction*/,
+         const WarmstartInformation& warmstart_information) {
+      this->compute_cauchy_step(subproblem, warmstart_information);
+      // TODO
+   }
+
+   // private member functions
+
+   void DoglegEvaluationSpace::compute_cauchy_step(const Subproblem& subproblem, const WarmstartInformation& warmstart_information) {
+      if (warmstart_information.objective_changed) {
+         // g^T g
+         this->objective_gradient_squared_norm = dot(this->objective_gradient, this->objective_gradient);
+         // B g = H_k ∇f(x_k)
+         subproblem.compute_hessian_vector_product(subproblem.current_iterate.primals.data(), this->objective_gradient.data(),
+            this->hessian_gradient_product.data());
+         // g^T B g = ∇f(x_k)^T H_k ∇f(x_k)
+         this->hessian_quadratic_product = dot(this->objective_gradient, this->hessian_gradient_product);
+         if (this->hessian_quadratic_product <= 0.) {
+            throw std::runtime_error("The objective Hessian is not positive definite");
+         }
+         // Cauchy step: d_C = - (g^T g)/(g^T B g) g
+         this->cauchy_step = this->objective_gradient;
+         const double scaling_factor = -this->objective_gradient_squared_norm / this->hessian_quadratic_product;
+         this->cauchy_step.scale(scaling_factor);
+      }
+   }
+
+   // find the positive real root to ax^2 + bx + c = 0
+   double DoglegEvaluationSpace::compute_positive_root_quadratic_equation(double a, double b, double c) {
+      // avoid catastrophic cancellation
+      const double delta = b*b - 4.*a*c;
+      if (delta < 0.) {
+         throw std::runtime_error("No real root");
+      }
+      // TODO check denominator
+      return (2.*c)/(-b - std::sqrt(delta));
+   }
+} // namespace

uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp

@@ -0,0 +1,51 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#ifndef UNO_DOGLEGEVALUATIONSPACE_H
+#define UNO_DOGLEGEVALUATIONSPACE_H
+
+#include "optimization/EvaluationSpace.hpp"
+#include "linear_algebra/Vector.hpp"
+#include "tools/Infinity.hpp"
+
+namespace uno {
+   // forward declaration
+   class Direction;
+   template <typename ElementType>
+   class SymmetricIndefiniteLinearSolver;
+
+   class DoglegEvaluationSpace: public EvaluationSpace {
+   public:
+      DoglegEvaluationSpace() = default;
+      ~DoglegEvaluationSpace() override = default;
+
+      // Newton step
+      Vector<double> objective_gradient{};
+      Vector<double> newton_step{};
+      double newton_step_squared_norm{INF<double>};
+      // Cauchy step
+      double objective_gradient_squared_norm{INF<double>};
+      Vector<double> hessian_gradient_product{};
+      double hessian_quadratic_product{INF<double>};
+      Vector<double> cauchy_step{};
+
+      void initialize_memory(const Subproblem& subproblem);
+
+      void evaluate_constraint_jacobian(const OptimizationProblem& /*problem*/, Iterate& /*iterate*/) override;
+      void compute_constraint_jacobian_vector_product(const Vector<double>& /*vector*/, Vector<double>& result) const override;
+      void compute_constraint_jacobian_transposed_vector_product(const Vector<double>& vector, Vector<double>& result) const override;
+      [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem,
+         const Vector<double>& vector) const override;
+
+      void evaluate_objective_gradient(const Subproblem& subproblem, const WarmstartInformation& warmstart_information);
+      void compute_newton_step(Statistics& statistics, const Subproblem& subproblem,
+         SymmetricIndefiniteLinearSolver<double>& linear_solver, Direction& direction, const WarmstartInformation& warmstart_information);
+      void compute_dogleg(const Subproblem& subproblem, Direction& direction, const WarmstartInformation& warmstart_information);
+
+   private:
+      void compute_cauchy_step(const Subproblem& subproblem, const WarmstartInformation& warmstart_information);
+      [[nodiscard]] static double compute_positive_root_quadratic_equation(double a, double b, double c);
+   };
+} // namespace
+
+#endif // UNO_DOGLEGEVALUATIONSPACE_H

uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp

@@ -0,0 +1,39 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#include "DoglegMethod.hpp"
+#include "ingredients/subproblem/Subproblem.hpp"
+#include "ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolverFactory.hpp"
+#include "options/Options.hpp"
+
+namespace uno {
+   DoglegMethod::DoglegMethod(const Options& options):
+         linear_solver_name(options.get_string("linear_solver")) {
+   }
+
+   void DoglegMethod::initialize_memory(const Subproblem& subproblem) {
+      this->linear_solver = SymmetricIndefiniteLinearSolverFactory::create(this->linear_solver_name);
+      this->linear_solver->initialize_hessian(subproblem);
+
+      this->evaluation_space.initialize_memory(subproblem);
+   }
+
+   void DoglegMethod::solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius,
+         const Vector<double>& /*initial_point*/, Direction& direction, const WarmstartInformation& warmstart_information) {
+      this->evaluation_space.evaluate_objective_gradient(subproblem, warmstart_information);
+
+      const double squared_trust_region_radius = std::pow(trust_region_radius, 2.);
+      // first try the Newton step. This is the solution if within the trust region
+      this->evaluation_space.compute_newton_step(statistics, subproblem, *this->linear_solver, direction, warmstart_information);
+      if (this->evaluation_space.newton_step_squared_norm <= squared_trust_region_radius) {
+         return;
+      }
+      // if the trust region constraint is violated, compute the dogleg path: the broken path between the Cauchy step
+      // and the Newton step
+      this->evaluation_space.compute_dogleg(subproblem, direction, warmstart_information);
+   }
+
+   [[nodiscard]] EvaluationSpace& DoglegMethod::get_evaluation_space() {
+      return this->evaluation_space;
+   }
+} // namespace

uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp

@@ -0,0 +1,36 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#ifndef UNO_DOGLEGMETHOD_H
+#define UNO_DOGLEGMETHOD_H
+
+#include <memory>
+#include <string>
+#include "../BoundConstrainedSolver.hpp"
+#include "DoglegEvaluationSpace.hpp"
+#include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
+
+namespace uno {
+   // forward declaration
+   class Options;
+
+   class DoglegMethod: public BoundConstrainedSolver {
+   public:
+      explicit DoglegMethod(const Options& options);
+      ~DoglegMethod() override = default;
+
+      void initialize_memory(const Subproblem& subproblem) override;
+
+      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector<double>& initial_point,
+         Direction& direction, const WarmstartInformation& warmstart_information) override;
+
+      [[nodiscard]] EvaluationSpace& get_evaluation_space() override;
+
+   protected:
+      const std::string& linear_solver_name;
+      std::unique_ptr<DirectSymmetricIndefiniteLinearSolver<double>> linear_solver;
+      DoglegEvaluationSpace evaluation_space{};
+   };
+} // namespace
+
+#endif // UNO_DOGLEGMETHOD_H