[WIP] Dogleg method

CMakeLists.txt

@@ -68,6 +68,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/constraint_relaxation_strategies/FeasibilityRestoration.cpp

@@ -179,7 +179,6 @@ namespace uno {
       subproblem_solver.solve(statistics, subproblem, trust_region_radius, this->initial_point, direction, current_evaluations,
          warmstart_information);
       ++this->number_subproblems_solved;
-      direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables()));
       this->initial_point.fill(0.);
       DEBUG3 << direction << '\n';
    }

uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp

@@ -76,7 +76,6 @@ namespace uno {
       this->initial_point.fill(0.);
       this->subproblem_solver->solve(statistics, subproblem, trust_region_radius, this->initial_point, direction,
          current_evaluations, warmstart_information);
-      direction.norm = norm_inf(view(direction.primals, 0, this->original_problem.get_number_original_variables()));
       DEBUG3 << direction << '\n';
       warmstart_information.no_changes();
    }

uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp

@@ -226,6 +226,7 @@ namespace uno {
          direction.primals[variable_index] = std::min(std::max(direction.primals[variable_index], this->lower_bounds[variable_index]),
             this->upper_bounds[variable_index]);
       }
+      direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables()));
       // gather the multipliers
       this->set_multipliers(subproblem.number_variables, direction.multipliers);
       LPSolver::compute_dual_displacements(subproblem, direction.multipliers);

uno/ingredients/subproblem_solvers/COOLinearSystem.cpp

@@ -3,7 +3,6 @@
 
 #include "COOLinearSystem.hpp"
 #include "ingredients/subproblem/Subproblem.hpp"
-#include "linear_algebra/Indexing.hpp"
 
 namespace uno {
    COOLinearSystem::COOLinearSystem(int solver_indexing): solver_indexing(solver_indexing) {

uno/ingredients/subproblem_solvers/EQPSolver.cpp

@@ -17,6 +17,8 @@ namespace uno {
          linear_solver(SymmetricIndefiniteLinearSolverFactory::create(options.get_string("linear_solver"))) {
    }
 
+   EQPSolver::~EQPSolver() = default;
+
    void EQPSolver::initialize_memory(const Subproblem& subproblem) {
       if (!subproblem.has_hessian_matrix()) {
          throw std::runtime_error("The subproblem does not have an explicit Hessian matrix and cannot be solved with a direct linear solver");
@@ -67,6 +69,7 @@ namespace uno {
       }
       // assemble the full primal-dual direction
       subproblem.assemble_primal_dual_direction(linear_system.solution, direction);
+      direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables()));
    }
 
    SolverWorkspace& EQPSolver::get_workspace() {

uno/ingredients/subproblem_solvers/EQPSolver.hpp

@@ -16,7 +16,7 @@ namespace uno {
    class EQPSolver: public SubproblemSolver {
    public:
       explicit EQPSolver(const Options& options);
-      ~EQPSolver() override = default;
+      ~EQPSolver() override;
 
       void initialize_memory(const Subproblem& subproblem) override;
 

uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp

@@ -107,6 +107,7 @@ namespace uno {
             direction.multipliers.upper_bounds[variable_index] = bound_multiplier;
          }
       }
+      direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables()));
       // gather the multipliers
       for (size_t constraint_index = 0; constraint_index < subproblem.number_constraints; constraint_index++) {
          direction.multipliers.constraints[constraint_index] = solution.row_dual[constraint_index];

uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp

@@ -12,6 +12,7 @@
 #include "QPSolver.hpp"
 #include "QPSolverFactory.hpp"
 #include "WoodburyEQPSolver.hpp"
+#include "dogleg/DoglegMethod.hpp"
 #include "ingredients/subproblem/Subproblem.hpp"
 #include "options/Options.hpp"
 #include "tools/Logger.hpp"
@@ -38,7 +39,16 @@ namespace uno {
       const Subproblem subproblem(problem, current_iterate, hessian_model, inertia_correction_strategy);
       // if no inequality constraint and no trust region, allocate EQP solver
       // temporary fix: this is set only in interior-point methods
+      /*if (!subproblem.has_inequality_constraints() && uses_trust_region && subproblem.is_hessian_positive_definite()) {
+         // use the dogleg method
+         DEBUG << "Trust-region and no inequality constraints in the subproblem, allocating a dogleg solver\n";
+         auto subproblem_solver = std::make_unique<DoglegMethod>(options);
+         subproblem_solver->initialize_memory(subproblem);
+         return subproblem_solver;
+      }
+      else*/
       if (!subproblem.has_inequality_constraints() && !uses_trust_region && options.get_string("inequality_handling_method") == "interior_point") {
+         // no trust region
          if constexpr (std::is_same_v<HessianType, LBFGSHessian>) {
             DEBUG << "No inequality constraints in the subproblem, allocating an EQP solver with L-BFGS Hessian\n";
             // the hessian_model we pass has type LBFGSHessian
@@ -54,7 +64,7 @@ namespace uno {
          }
       }
       // otherwise, allocate LP/QP solver, depending on the presence of curvature in the subproblem
-      if (!subproblem.has_curvature()) {
+      else if (!subproblem.has_curvature()) {
          if (subproblem.number_constraints == 0) {
             DEBUG << "No curvature and only bound constraints in the subproblem, allocating a box LP solver\n";
             auto subproblem_solver = std::make_unique<BoxLPSolver>();

uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp

@@ -74,6 +74,7 @@ namespace uno {
 
       // assemble the full primal-dual direction
       subproblem.assemble_primal_dual_direction(solution_diagonal_part, direction);
+      direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables()));
    }
 
    SolverWorkspace& WoodburyEQPSolver::get_workspace() {

uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp

@@ -0,0 +1,36 @@
+// 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 "optimization/Direction.hpp"
+#include "tools/Logger.hpp"
+
+namespace uno {
+   DoglegMethod::DoglegMethod(const Options& options):
+         workspace(options) {
+   }
+
+   void DoglegMethod::initialize_memory(const Subproblem& subproblem) {
+      this->workspace.initialize_memory(subproblem);
+   }
+
+   void DoglegMethod::solve(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius,
+         const Vector<double>& /*initial_point*/, Direction& direction, Evaluations& current_evaluations,
+         const WarmstartInformation& warmstart_information) {
+      // first try the Newton step. This is the solution if within the trust region
+      this->workspace.compute_newton_step(statistics, subproblem, direction, current_evaluations, warmstart_information);
+      if (this->workspace.newton_step_squared_norm <= std::pow(trust_region_radius, 2.)) {
+         DEBUG << "The Newton step is within the trust region\n";
+         direction.norm = norm_2(view(direction.primals, 0, subproblem.problem.get_number_original_variables()));
+         return;
+      }
+      // if the trust region constraint is violated, compute the dogleg path: the broken path between the Cauchy step
+      // and the Newton step
+      this->workspace.compute_dogleg(subproblem, trust_region_radius, direction, current_evaluations, warmstart_information);
+   }
+
+   [[nodiscard]] SolverWorkspace& DoglegMethod::get_workspace() {
+      return this->workspace;
+   }
+} // namespace

uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.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_DOGLEGMETHOD_H
+#define UNO_DOGLEGMETHOD_H
+
+#include "../SubproblemSolver.hpp"
+#include "DoglegWorkspace.hpp"
+
+namespace uno {
+   // forward declaration
+   class Options;
+
+   class DoglegMethod: public SubproblemSolver {
+   public:
+      explicit DoglegMethod(const Options& options);
+      ~DoglegMethod() override = default;
+
+      void initialize_memory(const Subproblem& subproblem) override;
+
+      void solve(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius,
+         const Vector<double>& initial_point, Direction& direction, Evaluations& current_evaluations,
+         const WarmstartInformation& warmstart_information) override;
+
+      [[nodiscard]] SolverWorkspace& get_workspace() override;
+
+   protected:
+      DoglegWorkspace workspace;
+   };
+} // namespace
+
+#endif // UNO_DOGLEGMETHOD_H

uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp

@@ -0,0 +1,116 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#include "DoglegWorkspace.hpp"
+#include "ingredients/subproblem/Subproblem.hpp"
+#include "optimization/Direction.hpp"
+#include "optimization/Iterate.hpp"
+#include "optimization/OptimizationProblem.hpp"
+#include "optimization/WarmstartInformation.hpp"
+#include "symbolic/Sum.hpp"
+#include "tools/Logger.hpp"
+
+namespace uno {
+   DoglegWorkspace::DoglegWorkspace(const Options& options):
+         newton_solver(options) {
+   }
+
+   void DoglegWorkspace::initialize_memory(const Subproblem& subproblem) {
+      // Newton solver
+      this->newton_solver.initialize_memory(subproblem);
+      this->initial_point.resize(subproblem.number_variables);
+      // Newton step
+      this->g.resize(subproblem.number_variables);
+      this->newton_step.resize(subproblem.number_variables);
+      // Cauchy step
+      this->Hg.resize(subproblem.number_variables);
+      this->cauchy_step.resize(subproblem.number_variables);
+   }
+
+   double DoglegWorkspace::compute_hessian_quadratic_form(const Subproblem& subproblem, const Vector<double>& vector) const {
+      return this->newton_solver.get_workspace().compute_hessian_quadratic_form(subproblem, vector);
+   }
+
+   void DoglegWorkspace::compute_newton_step(Statistics& statistics, const Subproblem& subproblem, Direction& direction,
+         Evaluations& current_evaluations, const WarmstartInformation& warmstart_information) {
+      if (warmstart_information.new_iterate) {
+         this->newton_solver.solve(statistics, subproblem, INF<double>, this->initial_point, direction, current_evaluations,
+            warmstart_information);
+         this->newton_step = direction.primals;
+         this->newton_step_squared_norm = dot(this->newton_step, this->newton_step);
+      }
+      DEBUG << "Newton step: " << this->newton_step << '\n';
+   }
+
+   void DoglegWorkspace::compute_dogleg(const Subproblem& subproblem, double trust_region_radius, Direction& direction,
+         Evaluations& current_evaluations, const WarmstartInformation& warmstart_information) {
+      this->compute_cauchy_step(subproblem, current_evaluations, warmstart_information);
+      DEBUG << "Cauchy step: " << this->cauchy_step << '\n';
+      double squared_norm_cauchy = dot(this->cauchy_step, this->cauchy_step);
+      if (trust_region_radius*trust_region_radius <= squared_norm_cauchy) {
+         DEBUG << "The Cauchy step is outside the trust region. Returning a clipped direction\n";
+         direction.primals = this->cauchy_step;
+         direction.primals.scale(trust_region_radius/std::sqrt(squared_norm_cauchy));
+         direction.norm = trust_region_radius;
+         DEBUG << "Scaled Cauchy direction: " << direction.primals << '\n';
+         return;
+      }
+
+      DEBUG << "Computing the dogleg step\n";
+      // define temporary vectors
+      Vector<double> u(subproblem.number_variables);
+      u = this->newton_step - this->cauchy_step;
+      Vector<double> v(subproblem.number_variables);
+      v = 2.*this->cauchy_step - this->newton_step;
+      const double a = dot(u, u);
+      const double b = 2.*dot(u, v);
+      const double c = dot(v, v) - trust_region_radius*trust_region_radius;
+      DEBUG << "(a, b, c) = " << a << ", " << b << ", " << c << '\n';
+      const double tau = DoglegWorkspace::compute_positive_root_quadratic_equation(a, b, c);
+      DEBUG << "tau = " << tau << '\n';
+      if (1. <= tau && tau <= 2.) {
+         direction.primals = this->cauchy_step + (1. - tau)*(this->newton_step - this->cauchy_step);
+         direction.norm = norm_2(view(direction.primals, 0, subproblem.problem.get_number_original_variables()));
+      }
+      else {
+         throw std::runtime_error("The dogleg step could not be computed");
+      }
+   }
+
+   // private member functions
+
+   void DoglegWorkspace::compute_cauchy_step(const Subproblem& subproblem, Evaluations& current_evaluations,
+         const WarmstartInformation& warmstart_information) {
+      if (warmstart_information.new_iterate) {
+         subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->g.data(), current_evaluations);
+         DEBUG << "g = " << this->g << '\n';
+         // gᵀ g
+         this->g_squared_norm = dot(this->g, this->g);
+         // H g
+         subproblem.compute_hessian_vector_product(subproblem.current_iterate.primals.data(), this->g.data(),
+            this->Hg.data());
+         DEBUG << "H g = " << this->Hg << '\n';
+         // gᵀ H g
+         this->gHg = dot(this->g, this->Hg);
+         DEBUG << "g^T H g = " << this->gHg << '\n';
+         if (this->gHg <= 0.) {
+            throw std::runtime_error("The Hessian is not positive definite");
+         }
+         // Cauchy step: d_C = - (gᵀ g)/(gᵀ B g) g
+         this->cauchy_step = this->g;
+         const double scaling_factor = -this->g_squared_norm / this->gHg;
+         this->cauchy_step.scale(scaling_factor);
+      }
+   }
+
+   // find the positive real root to ax^2 + bx + c = 0
+   double DoglegWorkspace::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/DoglegWorkspace.hpp

@@ -0,0 +1,54 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#ifndef UNO_DOGLEGWORKSPACE_H
+#define UNO_DOGLEGWORKSPACE_H
+
+#include "../SolverWorkspace.hpp"
+#include "ingredients/subproblem_solvers/EQPSolver.hpp"
+#include "linear_algebra/Vector.hpp"
+#include "tools/Infinity.hpp"
+
+namespace uno {
+   // forward declaration
+   class Direction;
+   class Evaluations;
+   class Options;
+   class Statistics;
+   class WarmstartInformation;
+
+   class DoglegWorkspace: public SolverWorkspace {
+   public:
+      explicit DoglegWorkspace(const Options& options);
+      ~DoglegWorkspace() override = default;
+
+      // Newton step
+      Vector<double> g{};
+      Vector<double> newton_step{};
+      double newton_step_squared_norm{INF<double>};
+      // Cauchy step
+      double g_squared_norm{INF<double>};
+      Vector<double> Hg{};
+      double gHg{INF<double>};
+      Vector<double> cauchy_step{};
+
+      void initialize_memory(const Subproblem& subproblem);
+
+      [[nodiscard]] double compute_hessian_quadratic_form(const Subproblem& subproblem, const Vector<double>& vector) const override;
+
+      void compute_newton_step(Statistics& statistics, const Subproblem& subproblem, Direction& direction,
+         Evaluations& current_evaluations, const WarmstartInformation& warmstart_information);
+      void compute_dogleg(const Subproblem& subproblem, double trust_region_radius, Direction& direction,
+         Evaluations& current_evaluations, const WarmstartInformation& warmstart_information);
+
+   private:
+      mutable EQPSolver newton_solver;
+      Vector<double> initial_point;
+
+      void compute_cauchy_step(const Subproblem& subproblem, Evaluations& current_evaluations,
+         const WarmstartInformation& warmstart_information);
+      [[nodiscard]] static double compute_positive_root_quadratic_equation(double a, double b, double c);
+   };
+} // namespace
+
+#endif // UNO_DOGLEGWORKSPACE_H