From 542e0df74c8c49eab59a9238158eb78d85b7ad72 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Sun, 16 Nov 2025 22:05:13 +0100 Subject: [PATCH 01/12] First draft of the structure of the dogleg method --- CMakeLists.txt | 1 + .../BoundConstrainedSolver.hpp | 32 +++++++ .../dogleg/DoglegEvaluationSpace.cpp | 88 +++++++++++++++++++ .../dogleg/DoglegEvaluationSpace.hpp | 50 +++++++++++ .../dogleg/DoglegMethod.cpp | 38 ++++++++ .../dogleg/DoglegMethod.hpp | 36 ++++++++ 6 files changed, 245 insertions(+) create mode 100644 uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp create mode 100644 uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp create mode 100644 uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp create mode 100644 uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp create mode 100644 uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index e05f4cab6..9d1cc0fe2 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp b/uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp new file mode 100644 index 000000000..d096d06cd --- /dev/null +++ b/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 + 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& initial_point, + Direction& direction, const WarmstartInformation& warmstart_information) override = 0; + + [[nodiscard]] virtual EvaluationSpace& get_evaluation_space() = 0; + }; +} // namespace + +#endif // UNO_TRUSTREGIONSOLVER_H \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp new file mode 100644 index 000000000..ccefecae6 --- /dev/null +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp @@ -0,0 +1,88 @@ +// 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/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_objective_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& /*vector*/, Vector& result) const { + result.fill(0.); + } + + void DoglegEvaluationSpace::compute_constraint_jacobian_transposed_vector_product(const Vector& /*vector*/, + Vector& result) const { + result.fill(0.); + } + + double DoglegEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, + const Vector& /*vector*/) const { + throw std::runtime_error("Not implemented yet"); + } + + void DoglegEvaluationSpace::compute_newton_step(const Subproblem& subproblem, SymmetricIndefiniteLinearSolver& /*linear_solver*/, + 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()); + // TODO compute Newton step + // linear_solver.solve_indefinite_system(statistics, subproblem, direction, warmstart_information); + 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_objective_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_objective_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 \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp new file mode 100644 index 000000000..d7f06a2a8 --- /dev/null +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp @@ -0,0 +1,50 @@ +// 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 + class SymmetricIndefiniteLinearSolver; + + class DoglegEvaluationSpace: public EvaluationSpace { + public: + DoglegEvaluationSpace() = default; + ~DoglegEvaluationSpace() override = default; + + // Newton step + Vector objective_gradient{}; + Vector newton_step{}; + double newton_step_squared_norm{INF}; + // Cauchy step + double objective_gradient_squared_norm{INF}; + Vector hessian_objective_product{}; + double hessian_quadratic_product{INF}; + Vector 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& /*vector*/, Vector& result) const override; + void compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const override; + [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, + const Vector& vector) const override; + + void compute_newton_step(const Subproblem& subproblem, SymmetricIndefiniteLinearSolver& linear_solver, + 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 \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp new file mode 100644 index 000000000..54d008864 --- /dev/null +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp @@ -0,0 +1,38 @@ +// 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& /*initial_point*/, Direction& direction, const WarmstartInformation& 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(subproblem, *this->linear_solver, warmstart_information); + if (this->evaluation_space.newton_step_squared_norm <= squared_trust_region_radius) { + // TODO save Newton step into direction + 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->linear_solver->get_evaluation_space(); + } +} // namespace \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp new file mode 100644 index 000000000..b7c2510a8 --- /dev/null +++ b/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 +#include +#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& 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> linear_solver; + DoglegEvaluationSpace evaluation_space{}; + }; +} // namespace + +#endif // UNO_DOGLEGMETHOD_H \ No newline at end of file From cabaad4b8b8afecc9e1ceba9cd2574f9c0fbd9ef Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 27 Nov 2025 15:12:59 +0100 Subject: [PATCH 02/12] Some progress. TODO: compute dogleg --- .../dogleg/DoglegEvaluationSpace.cpp | 25 +++++++++++-------- .../dogleg/DoglegEvaluationSpace.hpp | 7 +++--- .../dogleg/DoglegMethod.cpp | 9 ++++--- 3 files changed, 24 insertions(+), 17 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp index ccefecae6..c067f9a5f 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp @@ -4,6 +4,7 @@ #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" @@ -13,7 +14,7 @@ namespace uno { this->objective_gradient.resize(subproblem.number_variables); this->newton_step.resize(subproblem.number_variables); // Cauchy step - this->hessian_objective_product.resize(subproblem.number_variables); + this->hessian_gradient_product.resize(subproblem.number_variables); this->cauchy_step.resize(subproblem.number_variables); } @@ -35,13 +36,19 @@ namespace uno { throw std::runtime_error("Not implemented yet"); } - void DoglegEvaluationSpace::compute_newton_step(const Subproblem& subproblem, SymmetricIndefiniteLinearSolver& /*linear_solver*/, - const WarmstartInformation& warmstart_information) { + 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& 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()); - // TODO compute Newton step - // linear_solver.solve_indefinite_system(statistics, subproblem, direction, warmstart_information); + 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); } } @@ -60,17 +67,15 @@ namespace uno { 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_objective_product.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_objective_product); + 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; + const double scaling_factor = -this->objective_gradient_squared_norm / this->hessian_quadratic_product; this->cauchy_step.scale(scaling_factor); } } diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp index d7f06a2a8..da955a70c 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp @@ -25,7 +25,7 @@ namespace uno { double newton_step_squared_norm{INF}; // Cauchy step double objective_gradient_squared_norm{INF}; - Vector hessian_objective_product{}; + Vector hessian_gradient_product{}; double hessian_quadratic_product{INF}; Vector cauchy_step{}; @@ -37,8 +37,9 @@ namespace uno { [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const override; - void compute_newton_step(const Subproblem& subproblem, SymmetricIndefiniteLinearSolver& linear_solver, - const WarmstartInformation& warmstart_information); + void evaluate_objective_gradient(const Subproblem& subproblem, const WarmstartInformation& warmstart_information); + void compute_newton_step(Statistics& statistics, const Subproblem& subproblem, + SymmetricIndefiniteLinearSolver& linear_solver, Direction& direction, const WarmstartInformation& warmstart_information); void compute_dogleg(const Subproblem& subproblem, Direction& direction, const WarmstartInformation& warmstart_information); private: diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp index 54d008864..26e154407 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp @@ -18,13 +18,14 @@ namespace uno { this->evaluation_space.initialize_memory(subproblem); } - void DoglegMethod::solve(Statistics& /*statistics*/, Subproblem& subproblem, double trust_region_radius, + void DoglegMethod::solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector& /*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(subproblem, *this->linear_solver, warmstart_information); + 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) { - // TODO save Newton step into direction return; } // if the trust region constraint is violated, compute the dogleg path: the broken path between the Cauchy step @@ -33,6 +34,6 @@ namespace uno { } [[nodiscard]] EvaluationSpace& DoglegMethod::get_evaluation_space() { - return this->linear_solver->get_evaluation_space(); + return this->evaluation_space; } } // namespace \ No newline at end of file From 449390e05d5167a16017d4661bd3b359b2273e10 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 27 Feb 2026 13:27:47 +0100 Subject: [PATCH 03/12] Fixes --- .../BoundConstrainedSolver.hpp | 7 +-- .../dogleg/DoglegMethod.cpp | 12 ++--- .../dogleg/DoglegMethod.hpp | 11 ++--- ...valuationSpace.cpp => DoglegWorkspace.cpp} | 45 ++++++------------- ...valuationSpace.hpp => DoglegWorkspace.hpp} | 28 ++++++------ 5 files changed, 44 insertions(+), 59 deletions(-) rename uno/ingredients/subproblem_solvers/dogleg/{DoglegEvaluationSpace.cpp => DoglegWorkspace.cpp} (59%) rename uno/ingredients/subproblem_solvers/dogleg/{DoglegEvaluationSpace.hpp => DoglegWorkspace.hpp} (58%) diff --git a/uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp b/uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp index d096d06cd..cb2cdd505 100644 --- a/uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp +++ b/uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp @@ -22,10 +22,11 @@ namespace uno { void initialize_memory(const Subproblem& subproblem) override = 0; - void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector& initial_point, - Direction& direction, const WarmstartInformation& warmstart_information) override = 0; + void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, + const Vector& initial_point, Direction& direction, Evaluations& current_evaluations, + const WarmstartInformation& warmstart_information) override = 0; - [[nodiscard]] virtual EvaluationSpace& get_evaluation_space() = 0; + [[nodiscard]] virtual SolverWorkspace& get_evaluation_space() = 0; }; } // namespace diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp index 26e154407..0a34654ca 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp @@ -19,21 +19,21 @@ namespace uno { } void DoglegMethod::solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, - const Vector& /*initial_point*/, Direction& direction, const WarmstartInformation& warmstart_information) { - this->evaluation_space.evaluate_objective_gradient(subproblem, warmstart_information); - + const Vector& /*initial_point*/, Direction& direction, Evaluations& current_evaluations, + const WarmstartInformation& 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); + this->evaluation_space.compute_newton_step(statistics, subproblem, direction, *this->linear_solver, current_evaluations, + 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); + this->evaluation_space.compute_dogleg(subproblem, direction, current_evaluations, warmstart_information); } - [[nodiscard]] EvaluationSpace& DoglegMethod::get_evaluation_space() { + [[nodiscard]] SolverWorkspace& DoglegMethod::get_evaluation_space() { return this->evaluation_space; } } // namespace \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp index b7c2510a8..07b3087e8 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp @@ -7,7 +7,7 @@ #include #include #include "../BoundConstrainedSolver.hpp" -#include "DoglegEvaluationSpace.hpp" +#include "DoglegWorkspace.hpp" #include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp" namespace uno { @@ -21,15 +21,16 @@ namespace uno { void initialize_memory(const Subproblem& subproblem) override; - void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector& initial_point, - Direction& direction, const WarmstartInformation& warmstart_information) override; + void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, + const Vector& initial_point, Direction& direction, Evaluations& current_evaluations, + const WarmstartInformation& warmstart_information) override; - [[nodiscard]] EvaluationSpace& get_evaluation_space() override; + [[nodiscard]] SolverWorkspace& get_evaluation_space() override; protected: const std::string& linear_solver_name; std::unique_ptr> linear_solver; - DoglegEvaluationSpace evaluation_space{}; + DoglegWorkspace evaluation_space{}; }; } // namespace diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp similarity index 59% rename from uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp rename to uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp index c067f9a5f..caecebbc0 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp @@ -1,7 +1,7 @@ // Copyright (c) 2025 Charlie Vanaret // Licensed under the MIT license. See LICENSE file in the project directory for details. -#include "DoglegEvaluationSpace.hpp" +#include "DoglegWorkspace.hpp" #include "ingredients/subproblem/Subproblem.hpp" #include "ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp" #include "optimization/Direction.hpp" @@ -9,7 +9,7 @@ #include "optimization/WarmstartInformation.hpp" namespace uno { - void DoglegEvaluationSpace::initialize_memory(const Subproblem& subproblem) { + void DoglegWorkspace::initialize_memory(const Subproblem& subproblem) { // Newton step this->objective_gradient.resize(subproblem.number_variables); this->newton_step.resize(subproblem.number_variables); @@ -18,42 +18,25 @@ namespace uno { 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& /*vector*/, Vector& result) const { - result.fill(0.); - } - - void DoglegEvaluationSpace::compute_constraint_jacobian_transposed_vector_product(const Vector& /*vector*/, - Vector& result) const { - result.fill(0.); - } - - double DoglegEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, + double DoglegWorkspace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const Vector& /*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& linear_solver, Direction& direction, const WarmstartInformation& warmstart_information) { - if (warmstart_information.objective_changed) { + void DoglegWorkspace::compute_newton_step(Statistics& statistics, const Subproblem& subproblem, Direction& direction, + SymmetricIndefiniteLinearSolver& linear_solver, Evaluations& current_evaluations, + const WarmstartInformation& warmstart_information) { + if (warmstart_information.new_iterate) { // 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); + subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->objective_gradient.data(), + current_evaluations); + linear_solver.solve_indefinite_system(statistics, subproblem, direction, current_evaluations, 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*/, + void DoglegWorkspace::compute_dogleg(const Subproblem& subproblem, Direction& /*direction*/, Evaluations& /*current_evaluations*/, const WarmstartInformation& warmstart_information) { this->compute_cauchy_step(subproblem, warmstart_information); // TODO @@ -61,8 +44,8 @@ namespace uno { // private member functions - void DoglegEvaluationSpace::compute_cauchy_step(const Subproblem& subproblem, const WarmstartInformation& warmstart_information) { - if (warmstart_information.objective_changed) { + void DoglegWorkspace::compute_cauchy_step(const Subproblem& subproblem, const WarmstartInformation& warmstart_information) { + if (warmstart_information.new_iterate) { // g^T g this->objective_gradient_squared_norm = dot(this->objective_gradient, this->objective_gradient); // B g = H_k ∇f(x_k) @@ -81,7 +64,7 @@ namespace uno { } // find the positive real root to ax^2 + bx + c = 0 - double DoglegEvaluationSpace::compute_positive_root_quadratic_equation(double a, double b, double c) { + 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.) { diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp similarity index 58% rename from uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp rename to uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp index da955a70c..fd0c8fb14 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegEvaluationSpace.hpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp @@ -1,23 +1,26 @@ // 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 +#ifndef UNO_DOGLEGWORKSPACE_H +#define UNO_DOGLEGWORKSPACE_H -#include "optimization/EvaluationSpace.hpp" +#include "../SolverWorkspace.hpp" #include "linear_algebra/Vector.hpp" #include "tools/Infinity.hpp" namespace uno { // forward declaration class Direction; + class Evaluations; + class Statistics; template class SymmetricIndefiniteLinearSolver; + class WarmstartInformation; - class DoglegEvaluationSpace: public EvaluationSpace { + class DoglegWorkspace: public SolverWorkspace { public: - DoglegEvaluationSpace() = default; - ~DoglegEvaluationSpace() override = default; + DoglegWorkspace() = default; + ~DoglegWorkspace() override = default; // Newton step Vector objective_gradient{}; @@ -31,16 +34,13 @@ namespace uno { void initialize_memory(const Subproblem& subproblem); - void evaluate_constraint_jacobian(const OptimizationProblem& /*problem*/, Iterate& /*iterate*/) override; - void compute_constraint_jacobian_vector_product(const Vector& /*vector*/, Vector& result) const override; - void compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const override; [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector& vector) const override; - void evaluate_objective_gradient(const Subproblem& subproblem, const WarmstartInformation& warmstart_information); - void compute_newton_step(Statistics& statistics, const Subproblem& subproblem, - SymmetricIndefiniteLinearSolver& linear_solver, Direction& direction, const WarmstartInformation& warmstart_information); - void compute_dogleg(const Subproblem& subproblem, Direction& direction, const WarmstartInformation& warmstart_information); + void compute_newton_step(Statistics& statistics, const Subproblem& subproblem, Direction& direction, + SymmetricIndefiniteLinearSolver& linear_solver, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); + void compute_dogleg(const Subproblem& subproblem, Direction& direction, Evaluations& current_evaluations, + const WarmstartInformation& warmstart_information); private: void compute_cauchy_step(const Subproblem& subproblem, const WarmstartInformation& warmstart_information); @@ -48,4 +48,4 @@ namespace uno { }; } // namespace -#endif // UNO_DOGLEGEVALUATIONSPACE_H \ No newline at end of file +#endif // UNO_DOGLEGWORKSPACE_H \ No newline at end of file From fa7c4c6ac3bfd72dbea6acf0524f64d7214142f4 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 30 Mar 2026 13:38:12 +0200 Subject: [PATCH 04/12] Adapted to recent changes --- .../BoundConstrainedSolver.hpp | 33 ------------------- .../dogleg/DoglegMethod.cpp | 19 ++++++----- .../dogleg/DoglegMethod.hpp | 10 +++--- .../dogleg/DoglegWorkspace.cpp | 13 +++++--- .../dogleg/DoglegWorkspace.hpp | 7 ++-- 5 files changed, 28 insertions(+), 54 deletions(-) delete mode 100644 uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp diff --git a/uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp b/uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp deleted file mode 100644 index cb2cdd505..000000000 --- a/uno/ingredients/subproblem_solvers/BoundConstrainedSolver.hpp +++ /dev/null @@ -1,33 +0,0 @@ -// 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 - 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& initial_point, Direction& direction, Evaluations& current_evaluations, - const WarmstartInformation& warmstart_information) override = 0; - - [[nodiscard]] virtual SolverWorkspace& get_evaluation_space() = 0; - }; -} // namespace - -#endif // UNO_TRUSTREGIONSOLVER_H \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp index 0a34654ca..1c24b9c33 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp @@ -3,6 +3,7 @@ #include "DoglegMethod.hpp" #include "ingredients/subproblem/Subproblem.hpp" +#include "ingredients/subproblem_solvers/LinearSystem.hpp" #include "ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolverFactory.hpp" #include "options/Options.hpp" @@ -13,27 +14,29 @@ namespace uno { void DoglegMethod::initialize_memory(const Subproblem& subproblem) { this->linear_solver = SymmetricIndefiniteLinearSolverFactory::create(this->linear_solver_name); - this->linear_solver->initialize_hessian(subproblem); + auto& linear_system = this->linear_solver->get_linear_system(); + linear_system.initialize_hessian(subproblem); + this->linear_solver->initialize_memory(); - this->evaluation_space.initialize_memory(subproblem); + this->workspace.initialize_memory(subproblem); } - void DoglegMethod::solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, + void DoglegMethod::solve(Statistics& /*statistics*/, const Subproblem& subproblem, double trust_region_radius, const Vector& /*initial_point*/, Direction& direction, Evaluations& current_evaluations, const WarmstartInformation& 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, direction, *this->linear_solver, current_evaluations, + this->workspace.compute_newton_step(subproblem, direction, *this->linear_solver, current_evaluations, warmstart_information); - if (this->evaluation_space.newton_step_squared_norm <= squared_trust_region_radius) { + if (this->workspace.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, current_evaluations, warmstart_information); + this->workspace.compute_dogleg(subproblem, direction, current_evaluations, warmstart_information); } - [[nodiscard]] SolverWorkspace& DoglegMethod::get_evaluation_space() { - return this->evaluation_space; + [[nodiscard]] SolverWorkspace& DoglegMethod::get_workspace() { + return this->workspace; } } // namespace \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp index 07b3087e8..1feeb5094 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp @@ -6,7 +6,7 @@ #include #include -#include "../BoundConstrainedSolver.hpp" +#include "../SubproblemSolver.hpp" #include "DoglegWorkspace.hpp" #include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp" @@ -14,23 +14,23 @@ namespace uno { // forward declaration class Options; - class DoglegMethod: public BoundConstrainedSolver { + class DoglegMethod: public SubproblemSolver { 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, + void solve(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius, const Vector& initial_point, Direction& direction, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information) override; - [[nodiscard]] SolverWorkspace& get_evaluation_space() override; + [[nodiscard]] SolverWorkspace& get_workspace() override; protected: const std::string& linear_solver_name; std::unique_ptr> linear_solver; - DoglegWorkspace evaluation_space{}; + DoglegWorkspace workspace{}; }; } // namespace diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp index caecebbc0..93b0ab2c9 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp @@ -3,8 +3,10 @@ #include "DoglegWorkspace.hpp" #include "ingredients/subproblem/Subproblem.hpp" +#include "ingredients/subproblem_solvers/LinearSystem.hpp" #include "ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp" #include "optimization/Direction.hpp" +#include "optimization/Iterate.hpp" #include "optimization/OptimizationProblem.hpp" #include "optimization/WarmstartInformation.hpp" @@ -18,19 +20,22 @@ namespace uno { this->cauchy_step.resize(subproblem.number_variables); } - double DoglegWorkspace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, - const Vector& /*vector*/) const { + double DoglegWorkspace::compute_hessian_quadratic_form(const Subproblem& /*subproblem*/, const Vector& /*vector*/) const { throw std::runtime_error("Not implemented yet"); } - void DoglegWorkspace::compute_newton_step(Statistics& statistics, const Subproblem& subproblem, Direction& direction, + void DoglegWorkspace::compute_newton_step(const Subproblem& subproblem, Direction& direction, SymmetricIndefiniteLinearSolver& linear_solver, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information) { if (warmstart_information.new_iterate) { // g = ∇f(x_k) subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->objective_gradient.data(), current_evaluations); - linear_solver.solve_indefinite_system(statistics, subproblem, direction, current_evaluations, warmstart_information); + auto& linear_system = linear_solver.get_linear_system(); + linear_system.initialize_hessian(subproblem); + linear_solver.initialize_memory(); + linear_solver.solve_indefinite_system(linear_system.solution.data()); + subproblem.assemble_primal_dual_direction(linear_system.solution, direction); this->newton_step = direction.primals; this->newton_step_squared_norm = dot(this->newton_step, this->newton_step); } diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp index fd0c8fb14..3b801817e 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp @@ -34,11 +34,10 @@ namespace uno { void initialize_memory(const Subproblem& subproblem); - [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, - const Vector& vector) const override; + [[nodiscard]] double compute_hessian_quadratic_form(const Subproblem& subproblem, const Vector& vector) const override; - void compute_newton_step(Statistics& statistics, const Subproblem& subproblem, Direction& direction, - SymmetricIndefiniteLinearSolver& linear_solver, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); + void compute_newton_step(const Subproblem& subproblem, Direction& direction, SymmetricIndefiniteLinearSolver& linear_solver, + Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); void compute_dogleg(const Subproblem& subproblem, Direction& direction, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); From 7b36b585db4cbd54eab419397a92fe53635beb48 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 30 Mar 2026 14:19:14 +0200 Subject: [PATCH 05/12] Fixed subproblem solver factory --- .../subproblem_solvers/SubproblemSolverFactory.hpp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp index 8cb9c6611..9471e278b 100644 --- a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp +++ b/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,15 @@ 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 && options.get_string("inequality_handling_method") == "interior_point") { + 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(options); + subproblem_solver->initialize_memory(subproblem); + return subproblem_solver; + } + else if (!subproblem.has_inequality_constraints() && !uses_trust_region) { + // no trust region if constexpr (std::is_same_v) { 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 +63,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(); From f82d1b62d82b2f1746dea67d7af06f921b418132 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 30 Mar 2026 15:40:32 +0200 Subject: [PATCH 06/12] Form the Newton system --- .../subproblem_solvers/COOLinearSystem.cpp | 1 - .../dogleg/DoglegMethod.cpp | 1 + .../dogleg/DoglegWorkspace.cpp | 37 ++++++++++++------- .../dogleg/DoglegWorkspace.hpp | 6 +-- 4 files changed, 28 insertions(+), 17 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/COOLinearSystem.cpp b/uno/ingredients/subproblem_solvers/COOLinearSystem.cpp index 3a5398b6d..327c55a8e 100644 --- a/uno/ingredients/subproblem_solvers/COOLinearSystem.cpp +++ b/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) { diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp index 1c24b9c33..c34d9f943 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp @@ -17,6 +17,7 @@ namespace uno { auto& linear_system = this->linear_solver->get_linear_system(); linear_system.initialize_hessian(subproblem); this->linear_solver->initialize_memory(); + this->linear_solver->do_symbolic_analysis(); this->workspace.initialize_memory(subproblem); } diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp index 93b0ab2c9..4679c91ec 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp @@ -4,11 +4,12 @@ #include "DoglegWorkspace.hpp" #include "ingredients/subproblem/Subproblem.hpp" #include "ingredients/subproblem_solvers/LinearSystem.hpp" -#include "ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp" +#include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp" #include "optimization/Direction.hpp" #include "optimization/Iterate.hpp" #include "optimization/OptimizationProblem.hpp" #include "optimization/WarmstartInformation.hpp" +#include "tools/Statistics.hpp" namespace uno { void DoglegWorkspace::initialize_memory(const Subproblem& subproblem) { @@ -25,16 +26,26 @@ namespace uno { } void DoglegWorkspace::compute_newton_step(const Subproblem& subproblem, Direction& direction, - SymmetricIndefiniteLinearSolver& linear_solver, Evaluations& current_evaluations, + DirectSymmetricIndefiniteLinearSolver& linear_solver, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information) { if (warmstart_information.new_iterate) { - // g = ∇f(x_k) + // g = ∇f(xₖ) subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->objective_gradient.data(), current_evaluations); + // form matrix auto& linear_system = linear_solver.get_linear_system(); - linear_system.initialize_hessian(subproblem); - linear_solver.initialize_memory(); + Statistics statistics{}; + subproblem.evaluate_lagrangian_hessian(statistics, linear_system.matrix_values.data()); + std::cout << "MATRIX: " << linear_system.matrix_values << '\n'; + linear_solver.do_numerical_factorization(true); + const Inertia inertia = linear_solver.get_inertia(); + std::cout << "INERTIA = " << inertia << '\n'; + // form RHS + subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, linear_system.rhs.data(), current_evaluations); + linear_system.rhs.scale(-1.); + std::cout << "RHS: " << linear_system.rhs << '\n'; linear_solver.solve_indefinite_system(linear_system.solution.data()); + std::cout << "DOGLEG: NEWTON SOLUTION: " << linear_system.solution << '\n'; subproblem.assemble_primal_dual_direction(linear_system.solution, direction); this->newton_step = direction.primals; this->newton_step_squared_norm = dot(this->newton_step, this->newton_step); @@ -51,19 +62,19 @@ namespace uno { void DoglegWorkspace::compute_cauchy_step(const Subproblem& subproblem, const WarmstartInformation& warmstart_information) { if (warmstart_information.new_iterate) { - // g^T g + // gᵀ g this->objective_gradient_squared_norm = dot(this->objective_gradient, this->objective_gradient); - // B g = H_k ∇f(x_k) + // B g = Hₖ ∇f(xₖ) 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"); + // gᵀ B g = ∇f(xₖ)ᵀ Hₖ ∇f(xₖ) + this->hessian_directional_derivative = dot(this->objective_gradient, this->hessian_gradient_product); + if (this->hessian_directional_derivative <= 0.) { + throw std::runtime_error("The Hessian is not positive definite"); } - // Cauchy step: d_C = - (g^T g)/(g^T B g) g + // Cauchy step: d_C = - (gᵀ g)/(gᵀ B g) g this->cauchy_step = this->objective_gradient; - const double scaling_factor = -this->objective_gradient_squared_norm / this->hessian_quadratic_product; + const double scaling_factor = -this->objective_gradient_squared_norm / this->hessian_directional_derivative; this->cauchy_step.scale(scaling_factor); } } diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp index 3b801817e..605d2e506 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp @@ -14,7 +14,7 @@ namespace uno { class Evaluations; class Statistics; template - class SymmetricIndefiniteLinearSolver; + class DirectSymmetricIndefiniteLinearSolver; class WarmstartInformation; class DoglegWorkspace: public SolverWorkspace { @@ -29,14 +29,14 @@ namespace uno { // Cauchy step double objective_gradient_squared_norm{INF}; Vector hessian_gradient_product{}; - double hessian_quadratic_product{INF}; + double hessian_directional_derivative{INF}; Vector cauchy_step{}; void initialize_memory(const Subproblem& subproblem); [[nodiscard]] double compute_hessian_quadratic_form(const Subproblem& subproblem, const Vector& vector) const override; - void compute_newton_step(const Subproblem& subproblem, Direction& direction, SymmetricIndefiniteLinearSolver& linear_solver, + void compute_newton_step(const Subproblem& subproblem, Direction& direction, DirectSymmetricIndefiniteLinearSolver& linear_solver, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); void compute_dogleg(const Subproblem& subproblem, Direction& direction, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); From 2e5e623941abd2bc245d2b0331ee8deae2868c1c Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 30 Mar 2026 16:15:10 +0200 Subject: [PATCH 07/12] Embed an EQPSolver as Newton solver in DoglegWorkspace --- .../subproblem_solvers/EQPSolver.cpp | 2 + .../subproblem_solvers/EQPSolver.hpp | 2 +- .../dogleg/DoglegMethod.cpp | 16 ++------ .../dogleg/DoglegMethod.hpp | 7 +--- .../dogleg/DoglegWorkspace.cpp | 41 ++++++------------- .../dogleg/DoglegWorkspace.hpp | 10 +++-- 6 files changed, 26 insertions(+), 52 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/EQPSolver.cpp b/uno/ingredients/subproblem_solvers/EQPSolver.cpp index 7cdeb98bf..03d6de986 100644 --- a/uno/ingredients/subproblem_solvers/EQPSolver.cpp +++ b/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"); diff --git a/uno/ingredients/subproblem_solvers/EQPSolver.hpp b/uno/ingredients/subproblem_solvers/EQPSolver.hpp index bf621c47e..1f2013445 100644 --- a/uno/ingredients/subproblem_solvers/EQPSolver.hpp +++ b/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; diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp index c34d9f943..dfa55cbb5 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp @@ -3,32 +3,22 @@ #include "DoglegMethod.hpp" #include "ingredients/subproblem/Subproblem.hpp" -#include "ingredients/subproblem_solvers/LinearSystem.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")) { + workspace(options) { } void DoglegMethod::initialize_memory(const Subproblem& subproblem) { - this->linear_solver = SymmetricIndefiniteLinearSolverFactory::create(this->linear_solver_name); - auto& linear_system = this->linear_solver->get_linear_system(); - linear_system.initialize_hessian(subproblem); - this->linear_solver->initialize_memory(); - this->linear_solver->do_symbolic_analysis(); - this->workspace.initialize_memory(subproblem); } - void DoglegMethod::solve(Statistics& /*statistics*/, const Subproblem& subproblem, double trust_region_radius, + void DoglegMethod::solve(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius, const Vector& /*initial_point*/, Direction& direction, Evaluations& current_evaluations, const WarmstartInformation& 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->workspace.compute_newton_step(subproblem, direction, *this->linear_solver, current_evaluations, - warmstart_information); + this->workspace.compute_newton_step(statistics, subproblem, direction, current_evaluations, warmstart_information); if (this->workspace.newton_step_squared_norm <= squared_trust_region_radius) { return; } diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp index 1feeb5094..29e812120 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp @@ -4,11 +4,8 @@ #ifndef UNO_DOGLEGMETHOD_H #define UNO_DOGLEGMETHOD_H -#include -#include #include "../SubproblemSolver.hpp" #include "DoglegWorkspace.hpp" -#include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp" namespace uno { // forward declaration @@ -28,9 +25,7 @@ namespace uno { [[nodiscard]] SolverWorkspace& get_workspace() override; protected: - const std::string& linear_solver_name; - std::unique_ptr> linear_solver; - DoglegWorkspace workspace{}; + DoglegWorkspace workspace; }; } // namespace diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp index 4679c91ec..5f6d69afc 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp @@ -3,16 +3,19 @@ #include "DoglegWorkspace.hpp" #include "ingredients/subproblem/Subproblem.hpp" -#include "ingredients/subproblem_solvers/LinearSystem.hpp" -#include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp" #include "optimization/Direction.hpp" #include "optimization/Iterate.hpp" #include "optimization/OptimizationProblem.hpp" #include "optimization/WarmstartInformation.hpp" -#include "tools/Statistics.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); // Newton step this->objective_gradient.resize(subproblem.number_variables); this->newton_step.resize(subproblem.number_variables); @@ -25,31 +28,13 @@ namespace uno { throw std::runtime_error("Not implemented yet"); } - void DoglegWorkspace::compute_newton_step(const Subproblem& subproblem, Direction& direction, - DirectSymmetricIndefiniteLinearSolver& linear_solver, Evaluations& current_evaluations, - const WarmstartInformation& warmstart_information) { - if (warmstart_information.new_iterate) { - // g = ∇f(xₖ) - subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->objective_gradient.data(), - current_evaluations); - // form matrix - auto& linear_system = linear_solver.get_linear_system(); - Statistics statistics{}; - subproblem.evaluate_lagrangian_hessian(statistics, linear_system.matrix_values.data()); - std::cout << "MATRIX: " << linear_system.matrix_values << '\n'; - linear_solver.do_numerical_factorization(true); - const Inertia inertia = linear_solver.get_inertia(); - std::cout << "INERTIA = " << inertia << '\n'; - // form RHS - subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, linear_system.rhs.data(), current_evaluations); - linear_system.rhs.scale(-1.); - std::cout << "RHS: " << linear_system.rhs << '\n'; - linear_solver.solve_indefinite_system(linear_system.solution.data()); - std::cout << "DOGLEG: NEWTON SOLUTION: " << linear_system.solution << '\n'; - subproblem.assemble_primal_dual_direction(linear_system.solution, direction); - this->newton_step = direction.primals; - this->newton_step_squared_norm = dot(this->newton_step, this->newton_step); - } + void DoglegWorkspace::compute_newton_step(Statistics& statistics, const Subproblem& subproblem, Direction& direction, + Evaluations& current_evaluations, const WarmstartInformation& warmstart_information) { + Vector initial_point(0); + this->newton_solver.solve(statistics, subproblem, INF, initial_point, direction, current_evaluations, + warmstart_information); + this->newton_step = direction.primals; + this->newton_step_squared_norm = dot(this->newton_step, this->newton_step); } void DoglegWorkspace::compute_dogleg(const Subproblem& subproblem, Direction& /*direction*/, Evaluations& /*current_evaluations*/, diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp index 605d2e506..9fda17ad2 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp @@ -5,6 +5,7 @@ #define UNO_DOGLEGWORKSPACE_H #include "../SolverWorkspace.hpp" +#include "ingredients/subproblem_solvers/EQPSolver.hpp" #include "linear_algebra/Vector.hpp" #include "tools/Infinity.hpp" @@ -12,14 +13,13 @@ namespace uno { // forward declaration class Direction; class Evaluations; + class Options; class Statistics; - template - class DirectSymmetricIndefiniteLinearSolver; class WarmstartInformation; class DoglegWorkspace: public SolverWorkspace { public: - DoglegWorkspace() = default; + explicit DoglegWorkspace(const Options& options); ~DoglegWorkspace() override = default; // Newton step @@ -36,12 +36,14 @@ namespace uno { [[nodiscard]] double compute_hessian_quadratic_form(const Subproblem& subproblem, const Vector& vector) const override; - void compute_newton_step(const Subproblem& subproblem, Direction& direction, DirectSymmetricIndefiniteLinearSolver& linear_solver, + void compute_newton_step(Statistics& statistics, const Subproblem& subproblem, Direction& direction, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); void compute_dogleg(const Subproblem& subproblem, Direction& direction, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); private: + EQPSolver newton_solver; + 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); }; From 44887a098862020edab41cdaae000dc1861a5095 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 30 Mar 2026 16:21:34 +0200 Subject: [PATCH 08/12] Evaluate objective gradient (already evaluated for Newton step but stored in the linear system...) --- .../dogleg/DoglegWorkspace.cpp | 16 +++++++++++----- .../dogleg/DoglegWorkspace.hpp | 3 ++- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp index 5f6d69afc..aec4f0802 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp @@ -37,25 +37,31 @@ namespace uno { this->newton_step_squared_norm = dot(this->newton_step, this->newton_step); } - void DoglegWorkspace::compute_dogleg(const Subproblem& subproblem, Direction& /*direction*/, Evaluations& /*current_evaluations*/, + void DoglegWorkspace::compute_dogleg(const Subproblem& subproblem, Direction& /*direction*/, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information) { - this->compute_cauchy_step(subproblem, warmstart_information); + this->compute_cauchy_step(subproblem, current_evaluations, warmstart_information); // TODO } // private member functions - void DoglegWorkspace::compute_cauchy_step(const Subproblem& subproblem, const WarmstartInformation& warmstart_information) { + 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->objective_gradient.data(), + current_evaluations); + std::cout << "∇f(xₖ) = " << this->objective_gradient << '\n'; // gᵀ g this->objective_gradient_squared_norm = dot(this->objective_gradient, this->objective_gradient); // B g = Hₖ ∇f(xₖ) subproblem.compute_hessian_vector_product(subproblem.current_iterate.primals.data(), this->objective_gradient.data(), this->hessian_gradient_product.data()); + std::cout << "Hₖ ∇f(xₖ) = " << this->hessian_gradient_product << '\n'; // gᵀ B g = ∇f(xₖ)ᵀ Hₖ ∇f(xₖ) this->hessian_directional_derivative = dot(this->objective_gradient, this->hessian_gradient_product); - if (this->hessian_directional_derivative <= 0.) { - throw std::runtime_error("The Hessian is not positive definite"); + std::cout << "∇f(xₖ)ᵀ Hₖ ∇f(xₖ) = " << this->hessian_directional_derivative << '\n'; + if (this->hessian_directional_derivative < 0.) { + throw std::runtime_error("The Hessian is not positive semi-definite"); } // Cauchy step: d_C = - (gᵀ g)/(gᵀ B g) g this->cauchy_step = this->objective_gradient; diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp index 9fda17ad2..f93ea5cfd 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp @@ -44,7 +44,8 @@ namespace uno { private: EQPSolver newton_solver; - void compute_cauchy_step(const Subproblem& subproblem, const WarmstartInformation& warmstart_information); + 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 From 03ab2c385fb327f36391ba0b7ef3e7d2481d6406 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 30 Mar 2026 21:14:44 +0200 Subject: [PATCH 09/12] Working dogleg method + have the subproblem solvers compute the step norm in their own norm --- .../FeasibilityRestoration.cpp | 1 - .../NoRelaxation.cpp | 1 - .../hessian_models/ExactHessian.cpp | 2 +- .../subproblem_solvers/BQPD/BQPDSolver.cpp | 1 + .../subproblem_solvers/EQPSolver.cpp | 1 + .../subproblem_solvers/HiGHS/HiGHSSolver.cpp | 1 + .../subproblem_solvers/WoodburyEQPSolver.cpp | 1 + .../dogleg/DoglegMethod.cpp | 9 +- .../dogleg/DoglegWorkspace.cpp | 87 +++++++++++++------ .../dogleg/DoglegWorkspace.hpp | 15 ++-- 10 files changed, 79 insertions(+), 40 deletions(-) diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp index 09ef10000..5b45c08e9 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp +++ b/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'; } diff --git a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp index bb497c304..6b8a12738 100644 --- a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp +++ b/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(); } diff --git a/uno/ingredients/hessian_models/ExactHessian.cpp b/uno/ingredients/hessian_models/ExactHessian.cpp index 8d79f8862..9eba690e6 100644 --- a/uno/ingredients/hessian_models/ExactHessian.cpp +++ b/uno/ingredients/hessian_models/ExactHessian.cpp @@ -30,7 +30,7 @@ namespace uno { } bool ExactHessian::is_positive_definite() const { - return false; + return true;//false; } void ExactHessian::initialize_statistics(Statistics& /*statistics*/) const { diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp index ebe8a4443..ee2efe760 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp +++ b/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); diff --git a/uno/ingredients/subproblem_solvers/EQPSolver.cpp b/uno/ingredients/subproblem_solvers/EQPSolver.cpp index 03d6de986..3c3b129d5 100644 --- a/uno/ingredients/subproblem_solvers/EQPSolver.cpp +++ b/uno/ingredients/subproblem_solvers/EQPSolver.cpp @@ -67,6 +67,7 @@ namespace uno { direction.status = SubproblemStatus::INFEASIBLE; return; } + direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables())); // assemble the full primal-dual direction subproblem.assemble_primal_dual_direction(linear_system.solution, direction); } diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp index 34766ed31..dbfc5be4d 100644 --- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp +++ b/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]; diff --git a/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp b/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp index 639f2d638..979076d77 100644 --- a/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp +++ b/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp @@ -71,6 +71,7 @@ namespace uno { // compute the low-rank correction this->compute_low_rank_correction(subproblem, linear_system, solution_diagonal_part); + direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables())); // assemble the full primal-dual direction subproblem.assemble_primal_dual_direction(solution_diagonal_part, direction); diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp index dfa55cbb5..c8325cead 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp @@ -3,6 +3,8 @@ #include "DoglegMethod.hpp" #include "ingredients/subproblem/Subproblem.hpp" +#include "optimization/Direction.hpp" +#include "tools/Logger.hpp" namespace uno { DoglegMethod::DoglegMethod(const Options& options): @@ -16,15 +18,16 @@ namespace uno { void DoglegMethod::solve(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius, const Vector& /*initial_point*/, Direction& direction, Evaluations& current_evaluations, const WarmstartInformation& 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->workspace.compute_newton_step(statistics, subproblem, direction, current_evaluations, warmstart_information); - if (this->workspace.newton_step_squared_norm <= squared_trust_region_radius) { + 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, direction, current_evaluations, warmstart_information); + this->workspace.compute_dogleg(subproblem, trust_region_radius, direction, current_evaluations, warmstart_information); } [[nodiscard]] SolverWorkspace& DoglegMethod::get_workspace() { diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp index aec4f0802..4d649128b 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp @@ -7,6 +7,8 @@ #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): @@ -16,31 +18,63 @@ namespace uno { 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->objective_gradient.resize(subproblem.number_variables); + this->g.resize(subproblem.number_variables); this->newton_step.resize(subproblem.number_variables); // Cauchy step - this->hessian_gradient_product.resize(subproblem.number_variables); + this->Hg.resize(subproblem.number_variables); this->cauchy_step.resize(subproblem.number_variables); } - double DoglegWorkspace::compute_hessian_quadratic_form(const Subproblem& /*subproblem*/, const Vector& /*vector*/) const { - throw std::runtime_error("Not implemented yet"); + double DoglegWorkspace::compute_hessian_quadratic_form(const Subproblem& subproblem, const Vector& 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) { - Vector initial_point(0); - this->newton_solver.solve(statistics, subproblem, INF, initial_point, direction, current_evaluations, - warmstart_information); - this->newton_step = direction.primals; - this->newton_step_squared_norm = dot(this->newton_step, this->newton_step); + if (warmstart_information.new_iterate) { + this->newton_solver.solve(statistics, subproblem, INF, 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, Direction& /*direction*/, Evaluations& current_evaluations, - const WarmstartInformation& warmstart_information) { + 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); - // TODO + 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 u(subproblem.number_variables); + u = this->newton_step - this->cauchy_step; + Vector 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 @@ -48,24 +82,23 @@ namespace uno { 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->objective_gradient.data(), - current_evaluations); - std::cout << "∇f(xₖ) = " << this->objective_gradient << '\n'; + subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->g.data(), current_evaluations); + DEBUG << "g = " << this->g << '\n'; // gᵀ g - this->objective_gradient_squared_norm = dot(this->objective_gradient, this->objective_gradient); - // B g = Hₖ ∇f(xₖ) - subproblem.compute_hessian_vector_product(subproblem.current_iterate.primals.data(), this->objective_gradient.data(), - this->hessian_gradient_product.data()); - std::cout << "Hₖ ∇f(xₖ) = " << this->hessian_gradient_product << '\n'; - // gᵀ B g = ∇f(xₖ)ᵀ Hₖ ∇f(xₖ) - this->hessian_directional_derivative = dot(this->objective_gradient, this->hessian_gradient_product); - std::cout << "∇f(xₖ)ᵀ Hₖ ∇f(xₖ) = " << this->hessian_directional_derivative << '\n'; - if (this->hessian_directional_derivative < 0.) { - throw std::runtime_error("The Hessian is not positive semi-definite"); + 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->objective_gradient; - const double scaling_factor = -this->objective_gradient_squared_norm / this->hessian_directional_derivative; + this->cauchy_step = this->g; + const double scaling_factor = -this->g_squared_norm / this->gHg; this->cauchy_step.scale(scaling_factor); } } diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp index f93ea5cfd..e2f74f03d 100644 --- a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp @@ -23,13 +23,13 @@ namespace uno { ~DoglegWorkspace() override = default; // Newton step - Vector objective_gradient{}; + Vector g{}; Vector newton_step{}; double newton_step_squared_norm{INF}; // Cauchy step - double objective_gradient_squared_norm{INF}; - Vector hessian_gradient_product{}; - double hessian_directional_derivative{INF}; + double g_squared_norm{INF}; + Vector Hg{}; + double gHg{INF}; Vector cauchy_step{}; void initialize_memory(const Subproblem& subproblem); @@ -38,11 +38,12 @@ namespace uno { void compute_newton_step(Statistics& statistics, const Subproblem& subproblem, Direction& direction, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); - void compute_dogleg(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: - EQPSolver newton_solver; + mutable EQPSolver newton_solver; + Vector initial_point; void compute_cauchy_step(const Subproblem& subproblem, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); From efec9a687bab0ba9d32a7316bfbbf483eef462fc Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 30 Mar 2026 21:16:23 +0200 Subject: [PATCH 10/12] Disabled dogleg for now --- uno/ingredients/hessian_models/ExactHessian.cpp | 2 +- .../subproblem_solvers/SubproblemSolverFactory.hpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/uno/ingredients/hessian_models/ExactHessian.cpp b/uno/ingredients/hessian_models/ExactHessian.cpp index 9eba690e6..8d79f8862 100644 --- a/uno/ingredients/hessian_models/ExactHessian.cpp +++ b/uno/ingredients/hessian_models/ExactHessian.cpp @@ -30,7 +30,7 @@ namespace uno { } bool ExactHessian::is_positive_definite() const { - return true;//false; + return false; } void ExactHessian::initialize_statistics(Statistics& /*statistics*/) const { diff --git a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp index 9471e278b..e4a529e1b 100644 --- a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp +++ b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp @@ -39,14 +39,14 @@ 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()) { + /*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(options); subproblem_solver->initialize_memory(subproblem); return subproblem_solver; } - else if (!subproblem.has_inequality_constraints() && !uses_trust_region) { + else*/ if (!subproblem.has_inequality_constraints() && !uses_trust_region) { // no trust region if constexpr (std::is_same_v) { DEBUG << "No inequality constraints in the subproblem, allocating an EQP solver with L-BFGS Hessian\n"; From 103f2862a4c2a9469615c080b01e795cf371e2df Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 30 Mar 2026 22:06:41 +0200 Subject: [PATCH 11/12] Restored condition in SubproblemSolverFactory --- uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp index e4a529e1b..c17423925 100644 --- a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp +++ b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp @@ -46,7 +46,8 @@ namespace uno { subproblem_solver->initialize_memory(subproblem); return subproblem_solver; } - else*/ if (!subproblem.has_inequality_constraints() && !uses_trust_region) { + 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) { DEBUG << "No inequality constraints in the subproblem, allocating an EQP solver with L-BFGS Hessian\n"; From e734fa22b7f894971c2bbaffe42f50447fdeef16 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 30 Mar 2026 22:14:34 +0200 Subject: [PATCH 12/12] Fixed computation of direction norm --- uno/ingredients/subproblem_solvers/EQPSolver.cpp | 2 +- uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/EQPSolver.cpp b/uno/ingredients/subproblem_solvers/EQPSolver.cpp index 3c3b129d5..a510aafce 100644 --- a/uno/ingredients/subproblem_solvers/EQPSolver.cpp +++ b/uno/ingredients/subproblem_solvers/EQPSolver.cpp @@ -67,9 +67,9 @@ namespace uno { direction.status = SubproblemStatus::INFEASIBLE; return; } - direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables())); // 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() { diff --git a/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp b/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp index 979076d77..7c5f7d118 100644 --- a/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp +++ b/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp @@ -71,10 +71,10 @@ namespace uno { // compute the low-rank correction this->compute_low_rank_correction(subproblem, linear_system, solution_diagonal_part); - direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables())); // 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() {